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ABSTRACT 

A  number  of  computational  techniques  are  presented 
for  the  automatic  machine  calculation  of  the  structure 
and  properties  of  diatomic  molecules.   An  automatic 
procedure  is  described  which  uses  the  Boolean  operations 
of  a  binary  computer  for  constructing  N-electron  pure 
state  wave  functions  from  a  set  of  one-electron  functions. 
Use  of  such  automatic  procedures  permits  one  to  calculate, 
by  variational  methods,  electronic  wave  functions  from  a 
large  number  of  electronic  configurations  and  avoid  the 
use  and  limitations  of  the  self-consistent  field 
procedure.   The  problem  of  minimizing  a  function  of 
several  variables  was  Investigated  and  a  search  procedure 
which  fits  successive  quadratic  surfaces  to  calculated 
values  was  developed  and  programmed.   The  calculation  of 
tables  of  molecular  Integrals  of  Slater-type  orbltals, 
which  have  been  published  separately,  is  described.   The 
numerical  method  for  solving  the  radial  Schroedlnger 
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equation  for  bound  states  is  described  and  an  improved 
eigenvalue  correction  formula  is  derived.   An  analysis  of 
the  procedure  and  some  practical  means  for  improving  its 
convergence  and  accuracy  are  also  given.   The  technique 
was  programmed  as  a  computer  subroutine  and  used  for 
obtaining  nuclear  wave  functions  and  expectation  values 
of  properties  of  some  diatomic  molecules.   An  LCAO 
calculation  of  an  electronic  wave  function  for  H^  ,  with 
a  full  optimization  of  the  non-linear  parameters  by  the 
minimization  procedure  mentioned  above,  was  performed  for 
the  purpose  of  investigating  the  role  of  the  orbital 
exponents  and  the  applicability  of  the  Virial  Theorem, 
the  Hellman-Feymann  Theorem,  and  the  static  charge  cloud 
model.   The  results  lead  to  some  suggestions  which  will 
be  of  use  when  applying  similar  methods  to  many-electron 
systems.   The  vibrational-rotational  spectrum  of  Hp,   HD, 
and  Dp  was  calculated  with  two  electronic  wave  functions 
reported  by  other  workers  and  excellent  agreement  with 
experiment  was  obtained.   A  number  of  molecular  constants 
which  are  required  in  the  interpretation  of  radiofrequency 
spectra  are  also  given. 
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SOME  COMPUTATIONAL  METHODS  FOR  THE 
STUDY  OF  DIATOMIC  MOLECULES 

I   INTRODUCTION 

1.   Purpose 

Precise  measurements  of  observable  phenomena  resulting 
from  the  internal  structure  of  diatomic  molecules  can  be 
obtained  by  spectroscopic  methods.   Accurate  experimental 
techniques  have  been  developed,  over  a  long  period  of  time, 
for  obtaining  absorption  and  emission  spectra  from  the 
ultraviolet  to  the  far  infrared  regions.   The  interpreta- 
tion of  such  spectra  by  the  application  of  quantum- 
mechanical  principles,  with  simplifying  assumptions,  has 
led  to  well-established  hypotheses  concerning  the  particles 
involved  and  their  configuration  within  the  molecule. 
This  has  given  a  satisfactory  theoretical  explanation  of 
chemical  valence,  an  understanding  of  observed  physical 
and  chemical  properties  of  gases,  predictions  of  the 
existence  and  properties  of  molecules,  and  methods  for 
detecting  the  presence  of  molecules  and  rare  isotopic 
species,  particularly  in  the  planets,  comets,  and  upper 
atmosphere. 
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More  recently,  the  development,  by  N.  Ramsey  and 
others,  of  I.  I.  Rabl's  [l]  molecular  beam  magnetic  and 
electric  resonance  methods  has  given  extremely  precise 
measurements  in  the  low  energy  radiof requency  range  where 
the  spectrum  is  produced  by  relatively  weak  interactions 
between  the  nuclei,  the  surrounding  electronic  charge 
cloud,  and  external  magnetic  or  electric  fields. 

In  principle,  at  least,  all  of  these  observations 
can  be  predicted,  ab  initio,  from  hypotheses  concerning 
the  particles  involved  and  the  laws  of  quantum  mechanics 
but  the  mathematical  complexity  of  the  equations  for  even 
the  simplest  molecular  system  makes  their  solution  a 
formidable  task. 

The  relatively  recent  development  of  automatic  high- 
speed computers  has  made  it  possible  to  perform  some  of 
the  extensive  calculations  involved  and  there  is,  as  a 
result,  an  increasing  need  for  the  study,  development, 
and  application  of  appropriate  computational  methods. 
Among  the  problems  for  the  numerical  analyst,  one  finds 
many  frequently-encountered  computational  problems  such 
as  the  evaluation  of  rather  complicated  multi-dimensional 
definite  integrals,  matrix  eigenvalue  calculations,  the 
Integration  of  differential  equations,  and  minimum  value 
problems.   In  addition,  there  is  a  need  and  an  opportunity 
to  apply  non-numerical  procedures  to  the  analysis  of  the 


symmetry  of  molecules  and  the  mathematical  functions 
used  to  describe  their  structure.   Such  problems  can  be 
formulated  and  treated  by  techniques  which  use  the 
ability  of  the  electronic  computer  to  perform  logical 
operations. 

The  present  work  describes  some  computer-oriented 
techniques  for  the  application  of  the  laws  of  quantum 
mechanics  to  the  problems  of  determining  the  basic 
structure  of  diatomic  molecules.   At  times,  calculations 
of  the  type  described  here  yield  results  which  can  be 
compared  directly  with  experiment  to  check  the  validity 
of  the  hypotheses  and  techniques  used;  they  may  predict 
the  results  of  unobserved  experimental  phenomena;  and, 
in  many  Instances,  they  provide  quantities  which  are 
necessary  for  the  interpretation  of  experimental  data 
but  which  cannot  be  directly  measured. 


2.   The  Born-Oppenhelmer  Approximation 

A  calculation  of  the  wave  functions  of  a  molecule 
consisting  of  Interacting  electrons  and  nuclei  is 
simplified  by  the  use  of  the  Born-Oppenheimer  approxima- 
tion [2]  which  consists  in  treating  the  motion  of  the 
electrons  and  the  nuclei  separately.   Born  and 

Oppenheimer  have  shown  that  the  total  molecular  energy 

1/2 
can  be  expanded  in  powers  of   (m/M)  '       ,    where   m  is  the 

electronic  mass  and  M   is  an  average  nuclear  mass,  and 

that  separating  electronic  and  nuclear  motions  consists 

in  neglecting  high  order  terms  in  the  small  quantity 

1/2 
(m/M)  ^   .   Prom  another  point  of  view,  the  assumption 

is  that  the  nuclei,  being  so  heavy  as  compared  with  the 
electrons,  move  slowly  enough  for  the  electrons  to  adjust 
themselves  instantaneously  to  each  new  nuclear  configura- 
tion.  This  supports  the  further  assumption  that  the 
change  in  electron  configuration  during  nuclear  motion 
is  adiabatic:   i.e.,  there  is  a  negligible  probability  of 
a  change  from  one  electronic  state  to  another  during 
nuclear  motion. 

This  separation  of  the  various  motions  within  the 
molecule  separates  the  calculation  of  molecular  structure 
into  two  parts:   the  calculation  of  the  electronic  wave 
functions  as  described  in  Chapter  II,  and  the  relatively 
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simple  solution  of  the  one-dlmenslonal  radial  Schroedlnger 
equation  for  which  a  numerical  method  Is  given  In  Chapter 
IV. 

3.   Weak  Interactions 

In  addition  to  the  Born-Oppenhelmer  approximation, 
one  must  make  certain  other  simplifying  assumptions  to 
obtain  equations  which,  first  of  all,  actually  determine 
a  solution  and,  second,  lead  to  practical  means  of 
calculating  them.   Accordingly,  one  assumes  first  that 
the  nuclei  and  electrons  are  point  masses  and  charges 
obeying  the  Coulomb  law  of  force  and  that  the  magnetic 
moments  of  electron  spin,  nuclear  spin,  electron  orbital 
motion,  and  molecular  rotation  do  not  Interact.   In 
addition,  relatlvlstlc  effects  and  Interactions  with  an 
external  environment  are  omitted. 

These  effects  are  small  enough  to  allow  one  to  get 
useful  results  with  their  omission  and,  when  refinements 
are  needed,  to  use  first  and  second  order  perturbation 
methods  as  a  means  of  evaluating  the  quantities  omitted. 
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4 .   Atomic  Units 

Throughout  the  present  work,  all  quantities  will  be 
expressed  In  atomic  units  [5]  except  where  otherwise 
statedc   The  conversion  factors  used  are: 


1   a.u.  =  27.210  ev, 


=  -435  917  X10"-^°  erg 

=  219  474.62  cm"^ 

=  6.579  695x10^  Mc/sec. 

=  .52917  A 

=  274.078  Bohr  magnetons 
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II   THE  ELECTRONIC  WAVE  EQUATION 

1.   Introduction 

The  Hamlltonlan  of  the  energy  of  the  electrons  of  a 

s 

diatomic  molecule  with  the  nuclei   A  and  B  at  a  fixed 
internuclear  distance  R  ,  Is,  In  atomic  units. 


N       N-1    N 

(1.1)         '^  -     YZ\   +ZI  ZZ  ^i- 

1=1  ^        1=1  j  =  l+l    "^ 


where 


H.   =   -iA.  -Zr}-  Z,  r,  } 
1       2   1    a  ai    b  bi 


Here,   N   Is  the  number  of  electrons,   A.   Is  the  Laplaclan 

operator  for  the   1-th  electron,   Z   and   Z,   are  the 

charges  on  nuclei  A   and  B,   r  .   and   r,  ,   are  the 
'='  ai        bl 

distances  of  the   1-th  electron  from  A   and  B,  and   r. . 
Is  the  distance  between  the   1-th  and   J-th  electrons, 
A  wave  function  for  this  system  Is  an  N-electron 
function  JJ   ( 1,2,  .  .  ,  ,1,  .  ,,  .  ,N) ,  where   1   represents  the 
three  space  coordinates  x. ,  y. ,  z.,   and  the  spin 
coordinate   4.   of  the   1-th  electron.   It  Is  an  elgen- 
functlon  of  the  differential  equation 

(1.2)  %/  t  =   EI|J 
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for  some  eigenvalue   E  which  is  the  energy  of  the  system 
in  the  state  JJ  .   In  variational  methods,  approximate 
solutions  of  (1.2)  are  found  by  letting  a  trial  wave 
function  ^  range  over  some  necessarily  limited  set  of 
functions  and  taking  as  the  "best"  solutions  those  for 
which  the  expectation  value  of  the  energy, 

(1.5)     E  =     J    IJJ*%^  i  dV  /  r  ^*  JJ  dV  , 

attains  a  minimum.   Here,   /  dV  denotes  integration  over 
all   3N   space  and  N  spin  coordinates. 

An  important  aspect  of  the  computational  problem  is 
the  selection  of  trial  wave  functions  suitable  for  the 
problem  and  for  the  particular  eigenfunction  one  wants  to 
approximate. 

2.   Historical  Survey 

Hartree ' s  procedure  [4],  in  obtaining  solutions  of 
this  type  of  problem  for  atoms  was  to  form  a  simple 
product  trial  wave  function 


2.1)  JJ  =   ^^(1)  f^{2)    .     .     .    ^j^(N) 


where 


^l     =   ^^a  ,  ij/^     =      (j)^P  ,  tp^     =   (fj^a  ,  ip^     =   ^2^  , 
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with  a     and  P  denoting  the  one-electron  spin  elgen- 
functlons.   The  space  functions  were.  In  spherical  polar 
coordinates, 

(2.2)  ^.   =   R.(r)  Y,_^^_(e,^) 

where  the  Y.      's  are  the  complex  spherical  harmonics 
i! .  ,m. 

and  the   R. (r)  's  are  numerical  radial  functions.   By  a 
somewhat  heuristic  argument,  he  determined  the  radial 
functions  R. (r)   by  considering  the  Hamiltonian  of  the 
1-th  electron  in  the  potential  field  of  all  the  others. 


{2,3)        ^^      =      -  A  A^  -  Z  r-^  +  y~  f  y/J (  J  )  i^*{  j  )    r-^  dv^ 


where   Z   is  the  nuclear  charge  and   /  dv .   denotes 

u  J 

integration  over  the  space  and  spin  coordinates  of 
electron   j  . 

Slater  pointed  out  [5]  that  a  trial  wave  function 
should  satisfy  the  Paull  exclusion  principle  and  should 
be  antisymmetric  in  the  coordinates  of  the  electrons.   To 
accomplish  this,  he  suggested  the  use  of  "Slater 
Determinants"  ( SD ' s ) , 
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t  =   (n:)'^/^  det  [^^(1),.  .  .,^j^(N)] 


2A 


=       (NJ) 


-1/2 


^^(1)  ...  ^^(1) 


^^(N)  ...  ^j^(N) 


Fock  [6]  then  gave  a  rigorous  derivation  of  the  conditions 
under  which  the  energy  (1.3)  attains  a  minimum  when  computed 
with  the  antisymmetrlzed  product  wave  function  (2,4).   His 
result  was  that  the  one-electron  functions  which  diagonalize 
the  one-electron  operator, 


(2.5)  ^^    = 


-1 


N 


f  A.-Z  r.-^l 


J=l 


^^(j)  rT^d-P   )  ^  (j)  dv 


will  minimize  (l.^)-   When  (^ .      operates  on  a  one-electron 
function  i//.(i),  the  operator   P.  .   exchanges  the 
coordinates  of  electrons   i   and   j   before  the  integration 
takes  place.   Thus,  the  use  of  the  antisymmetrlzed  trial 
wave  functions  introduces  an  "exchange  term"  in  the  one- 
electron  operator. 

Integrating  over  angular  and  spin  coordinates  gives 
integro-dif f erential  equations  in  the  radial  functions 
R. (r)   corresponding  to  (2.3)  or  (2.5).   The  computational 
procedure  consists  in  forming  a  first  estimate  of  all 
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R. (r)  ' s  to  obtain  a  one-electron  effective  potential 
for  each  electron.   The  resulting  one-dlmenslonal 
differential  equation  Is  then  solved  for  a  new  set  of 
R. (r).   This  process  Is  repeated  until  the  Input   R. (r)  's 
are  the  same  as  the  output.   In  general,  the  use  of  such 
an  Iterative  procedure  with  a  single-determinant  wave 
function  Is  known  as  the  Hartree-Fock  (HF)  scheme  or  the 
"self-consistent  field"  ( SCF )  method. 

The  above  methods  are  characterized  by  the  fact 
that  an  N-electron  wave  function  Is  constructed  of  one- 
electron  functions.   This  Is  somewhat  distinct  from  other 
rather  successful  methods  In  which  the  N-electron  wave 
function  Is  formed  from  functions  of  the  coordinates  of 
several  electrons.   An  example  Is  the  treatment  of  the 
two-electron  He-like  systems  by  Hylleraas  .   Combinations 
of  the  two  points  of  view.  In  which  one-electron 
functions  are  combined  with  two-electron  functions,  have 
also  been  applied  [8].   The  great  success  of  Hylleraas 
on  two-electron  systems  tended  to  discourage  somewhat  the 
use  of  the  one-electron  approach  but  the  lack  of  ready 
applicability  of  the  Hylleraas-type  treatment  to  many- 
electron  systems  contributed  to  a  reviving  Interest  In  the 


For  a  brief  summary  and  bibliography  of  Hylleraas 
work,  see  Ref.  [TI* 


-  17  - 


latter.   Of  perhaps  more  Importance  in  this  revival  is 
the  refinement  of  mathematical  and  computational 
techniques  and  the  availability  of  electronic  computers 
capable  of  implementing  them. 

In  going  from  atoms  to  molecules,  the  Heitler- 
London-Pauling-Slater  or  valence  bond  (VB)  method  [7] 
was  applied  with  success  in  giving  a  quantum-mechanical 
picture  of  chemical  valence.   By  this  method,  the 
N-electron  molecular  wave  functions  are  constructed  in  a 
rather  simple  manner  from  the  wave  functions  of  the 
individual  atoms  in  the  molecule.   An  approach  to  the 
problem  from  another  point  of  view  is  contained  in  the 
Hund-Mulliken  or  molecular  orbital  (MO)  method,  which 
assigns  each  electron  to  a  one-electron  molecular  wave 
function  in  a  determinantal  molecular  wave  function. 
The  latter  has  its  advantages  in  describing  processes  of 
ionization  and  excitation.   The  line  of  development 
represented  by  the  very  successful  James  and  Coolidge  [9] 
treatment  of  the  hydrogen  molecule  with  a  two-electron 
function  falls  in  the  same  category  here  as  the  Hylleraas 
treatment  does  in  the  case  of  atoms. 

The  main  extension  of  the  Hartree-Fock  SCF  method  to 
molecules  is  in  the  treatment  given  by  Lennard-Jones 
[10,  11],  Hall  [12],  Roothaan  [I5],  and  others.   Here, 
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Instead  of  treating  numerical  functions,  a  convenient 
set  of  analytic  basis  functions,  "7-,,  7(p*  •  •  •  ?  °^  one- 
electron  space  coordinates  is  used  to  form  MO ' s  of  the 
form 


(2.6)  ^^  =  jz:  \c^,i 


Multiplying  all  MO ' s  by  a   and   P   spin  functions  gives 
twice  as  many  molecular  spin  orbitals  (MSO's).   These 
will  be  denoted  by 


•(//.   =   6^      a  if   i   is  odd 

%i  i+1 ) 
2.7)  2 


^1   =   <t>i   P 


if   i   is  even 


i^ 


By  forming  an  N-electron  SD  wave  function  of  the  form 
(2.4)  from  the  functions  (2.7),  the  integro-dif f erential 
equations  of  the  SCF  method  are  replaced  by  vector-matrix 
equations  which  determine  the  coefficients  in  (2.6).   A 
large  number  of  calculations  of  the  ground  states  of 
closed-shell  systems  have  been  performed  by  these  methods 
using,  as  basis  functions,  the  hydrogenic-like  atomic 
functions  to  be  defined  later. 

However,  the  HF  scheme  is  somewhat  difficult  or 
unsatisfactory  when  applied  to  open-shell  systems  [l4,  I5] 
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Furthermore,  even  in  the  latter  applications,  since  the 
spin-paired  MSO ' s  are  constrained  in  having  the  same 
space  functions,  the  motions  of  the  electrons  occupying 
them  are  not  correlated  by  the  Coulomb  repulsion  between 
them.   The  energy  of  this  constraint  has  come  to  be  known 
as  the  "correlation  energy"  and  has  been  calculated  to 
amount  to  approximately  1.0  to  1.6  ev.  per  electron  pair 
in  He-like  atoms  [l6].   Some  procedures,  such  as  the 
extended  Hartree-Fock  methods  [17],  have  given  some 
improvement  in  results. 

Another  difficulty  in  the  use  of  the  single-determinant 
trial  wave  function  is  that  it  is  suitable  only  for  certain 
eigenfunctions  of  the  problem  and  not  at  all  for  others. 
This  aspect  will  be  discussed  in  greater  detail  below.   The 
recent  line  of  development  which  has  given  some  improvement 
in  energy  values  and  properties  is  the  use  of  configuration 
interaction  (CI).   This  receives  some  mathematical 
encouragement  from  the  fact  that  the  set  of  SD ' s  obtained 
by  selecting  N   one-particle  functions  in  all  possible 
ways  from  a  complete  set  of  such  functions  is  complete  in 
the  Hilbert  space  of  all  antisymmetric   N-particle  functions. 
This  leads  one  to  believe  that  if  a  finite  set  of  one- 
particle  functions  is  sufficiently  adequate  in  describing 
the  motions  of  the  individual  electrons  in  a  system,  then 
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a  set  of  SD ' s  formed  from  them  should  be  able  to  represent 

the  N-electron  wave  functions  of  the  system.   One  rather 

successful  procedure  for  doing  this  is  to  go  through  the 

SCF  process  to  optimize  a  single-determinant  wave  function 

l|j   by  determining  the  "best"  MO '  s  of  the  form  (2.6).   In 
~o 

diagonalizing  the  one-electron  operator,  one  obtains 

extra  solutions,  or  "virtual  orbitals"  which  do  not  appear 

in  the  determinant.   Then,  the  orbital  energies,  which 

appear  as  eigenvalues  of  the  one-electron  operator,  may 

be  used  as  criteria  for  replacing  occupied  MO ' s  by  virtual 

MO '  s  in  various  ways.   The  set  of  SD's  l|j  ,  l|j  ,  .  .  .  ,  so 

-1  -2 

formed,  is  then  used  in  a  wave  function  of  the  form 

(2.8)  i  =  rZ  C;,  III  . 

K=0   ^   K 

By  the  variation  theorem,  the  optimal   C„  ' s  are  given  by 
the  solutions  of  the  secular  equations 


(2.9)    JZZ   (Hj^L  "  ^  ^KL^  ^L  =  °  ^        ^  =  0,1,2, 


L 


where 


"kl  =  J      *K^*L^^ 


* 


\l    =    J    \   \  ^v 
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In  this  way,  most  of  the  solution  Is  given  by  a  single 
determinant  and  the  mlxlng-ln  of  the  SD ' s  containing 
virtual  orbltals  gives  a  small  correction.   In  calcula- 
tions for  HF   Karo  and  Allen  [l8]  got  a  binding  energy 
of   1.3699  ev.   by  the  SCF  procedure  and  raised  this  to 
2.1166  ev.   by  using  seven  configurations.   The 
experimental  binding  energy  of   HF   is   6.O8  ev.   Kotani 
et  al.  [19]  raised  the  computed  binding  energy  of   0^ 
from   .795  ev.  to   3.629  ev.   by  using  I5  configurations. 
Here,  the  experimental  value  is   5.O8  ev. 

The  purpose  of  the  present  investigation  is  to  study 
the  possibility  of  directly  forming  trial  wave  functions 
of  the  form  (2,8)  from  MO ' s  which  are  not  necessarily 
SCF  MO's.   By  taking  all  possible  determinants  in  linear 
combinations  suitable  for  the  particular  eigenfunction 
desired,  one  does  not  impose  the  constraints  leading  to 
the  correlation  error  and  one  can  form  good  approximate 
wave  functions  for  eigenfunctlons  other  than  the  ground 
state. 

3.   Pure  State  Wave  Functions 

Considered  as  an  operator  on  the  Hilbert  space  of 
all  acceptable  wave  functions,  the  Hamiltonian  >/  , 
defined  by  (1.1),  is  a  Hermltlan  operator.   Here,  trial 
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wave  functions  for  the  variational  procedure  for 
calculating  approximations  to  the  eigenfunctions  and 
eigenvalues  of  y^    are  constructed  so  that  they  are, 
for  all  values  of  the  variational  parameters,  eigen- 
functions of  certain  Hermitian  operators  which  commute 
with  ')4     and  with  each  other.   A  wave  function  which  is 
simultaneously  an  eigenfunction  of  all  the  commuting 
operators  described  below,  but  not  necessarily  of  ^  , 
will  be  called  a  "pure  state  wave  function"  and  will  be 
classified  by  giving  the  corresponding  eigenvalues  of 
the  commuting  operators.   The  properties  of  these 
operators  imply  that  the  eigenfunctions  of  yr    corres- 
ponding to  each  eigenvalue  E  of  y^t    are  linear 
combinations  of  pure  state  wave  functions.   It  will  now 
be  shown  how  the  special  properties  of  the  pure  state 
wave  functions  may  be  used  to  advantage  in  the 
computational  procedure. 

Suppose  9^  is  one  of  the  commuting  operators  and 
that   A  and   A'   are  distinct  eigenvalues  of  f/l   .      Then, 
the  corresponding  eigenfunctions,  ^  and  ^  ,  of  ^ , 
whether  they  be  exact  eigenfunctions  of  Jr    or  simply 
variational  trial  wave  functions,  satisfy 


(3.1)  r  J""  _^'  dV 


0 
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and 

(3.2)  r  J  5V  $'  dv  =  0 , 

If  a  finite  set  of  eigenfunctions,  ^  ,  "^  ,  .  .  ,  ,  "^   of    09l 

-1  -2      -M       ^ 

is  then  chosen  and  a  trial  wave  function  is  expressed  in 
the  form 

(3.3)  $  =  y~  c  ^  , 

K    ^   K 
the  optimal  coefficients  of  (3.3)  satisfy 


(3.4)  YZ    (Hj,L  -  S^L  E)  ^L   =   0 

Li 

where 


■X- 


(3-5)  %L  ==   /  5  ^  i  dV 

^    K     L/ 


.* 


(3.6)  S^.   =   /   (1)   ({)  dV 


KL 


■-'  -Y. 


From  (3.1)  and  (3.2),  it  is  seen  that  the  J   's  can  be 
ordered  so  that  the  matrix  of  coefficients  of  the   C^  's 
in  (3.4)  appears  as  a  block  diagonal  matrix  with  each  block 
corresponding  to  a  distinct  eigenvalue  of  /^  .   Thus, 
each  eigenfunction  is  obtained  as  an  eigenvector  of  a 
particular  block  matrix  with  the  sum  in  (3-3)  going  over 
only  eigenfunctions  of  7/1     having  the  eigenvalue 


-  24  - 


corresponding  to  the  block  In  question.   In  this  way,  an 
approximation  to  the  exact  elgenfunctlon  for  the  selected 
value  of   A   is  found  by  solving  a  relatively  small  matrix 
equation.   Some  further  advantages  are  that:   the 
degeneracies  in  ^  associated  with  '9^  result  in  blocks 
on  the  diagonal  which  are  identical  so  that,  in  effect,  a 
number  of  different   A-state  solutions  may  be  obtained  in 
a  single  process;  non-linear  parameters  in  the  trial  wave 
functions  may  be  varied  without  destroying  the  properties 
(3"l)  and  (3-2);  and  the  variation  theorem  implies  that 
the  energies  obtained  are  upper  bounds  of  the  exact 
energies  for  the  particular   A-state  under  consideration. 

We  will  now  consider  specific  operators  which 
commute  with  /4  .      Denoting  the  molecular  axis  by   z,  the 
z-component  of  total  angular  momentum  is 

x=l 

where   i  .   =  -  y"^  S/StJ).   is  the  operator  corresponding 
to  the   z-component  of  angular  momentum  of  the   1-th 
electron.   Eigenfunctions  J(L  )   of  (^        corresponding 
to  the  eigenvalues  L   =0,  ±1,  ±2,.,,,  will  be  denoted 

by 

(3.8)  EI  ^    TT"  ,   Zf  ,..»  , 
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respecltvely,  in  spectroscopic  notation.   Following 
international  nomenclature,  we  let 


3.9)  A-  |L 


z 


The  second  operator  to  be  considered  is  (J^  ,    the 
reflection  of  the  wave  function  in  a  plane  containing 
the  molecular  axis.   It  can  easily  be  shown  from  the 
symmetry  of  the  molecule  that  the  dependence  of  l|J(L„) 
on  (j).  ,   i  =  1,2,...,N   is  such  that 

(3.10)         O-^l^L^)   =   1*L^)   =   i(-L2^  ' 

Therefore,  when  Jy      ^      0,  only  the  eigenfunctions 
I^(yV)   (i.e.  with  L   >   0)  need  be  found  since,  if 
needed,   _$_( -_/\,)   may  be  formed  by  simply  operating  with 
ff  .      It  should  be  noted  here  that  some  workers  use  the 
"  ±  "  in  the  notation  of  (3.8)  to  denote  the  two  eigen- 
functions 

(3.11)  ^"(A.)    =    lJ(y\.)  ±  t^.A.) 


of  cr   .      From  (3. 10)  and  (3,11)  it  is  seen  that  the 
eigenvalues  of  ty       are  ±1   and  that  the  eigenfunctions 
of  rf       are  the  real  and  imaginary  parts  of  the  eigen- 
functions (3-8)  of   .C  .   The  ^    -type  functions  are 
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doubly  degenerate  eigenfunctlons  of   ^C   and,  in  this 
case,  the  eigenf unctions  (3. 11)  of  CT^,  with  _/\_   =  0, 


will  be  constructed  and  denoted  by  ^ 


.+ 


In  the  special  case  of  a  homopolar  diatomic  molecule 
(having  nuclei  with  equal  charges),  another  operator 

r 

resulting  from  the  symmetry  of  the  molecule  is  -^  ,  the 
inversion  of  the  molecule  through  its  center.   This  can 
do  no  more  than  change  the  sign  of  the  wave  function,  so 
its  eigenvalues  are   ±1.   The  corresponding  eigenfunctlons 
are  referred  to  as  "gera^e"  and  "ungerade"  functions, 
respectively,  and  are  given  the  subscripts   g   and  u. 
Finally,  the  total  spin  operator  is  defined  by 

N 


S^^ 


•;3»i2)  X  =   HI  -^1 

1=1     ^ 

where  ^.       is  the  spin  operator  for  the   i-th  electron 
Since  ff    and  the  operators  considered  above  are 
independent  of  electron  spin,  they  commute  with  both 
and  the  component  ^„  x  ,    of  ^   along  any  axis   z'o   The 
eigenvalues  of  jt)         are   s(s  +  l),  where 

s   =   0,  1,  ...,  N/2  if   N   is  even 

(3.13) 

s   =   1/2,  3/2,..,,  N/2       if   N   is  odd  , 

ard  are  of  degeneracy   2s +1.   This  degeneracy  is  rem.c\  _ 
t'Y   J^,  i   which  has  the  eigenvalues 
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(3.1^)         s^,  =   -s,  -s+1,...,  s-1,  s   . 

If  yV  "f      0,  one  usually  considers  electronic  wave 

functions  whose  electron  spin  is  coupled  to  the  angular 

momentum  oL     •      In  this  case,  the   z  -axis  is  taken  as 
z 

the  molecular  axis,  z,  and   s    is  denoted,  in  the 


standard  notation  of  the  literature,  by  \ (the 

quantum  number  )     is  not  to  be  confused  with  the  ^ 


eigenfunction  in  (3o8)),   Therefore,  the  prime  will  be 
dropped  from  j<3  , ,  but  one  is  to  keep  in  mind  that  the 
spin  axis  is  arbitrary  and  other  coupling  schemes  may 
be  considered. 

In  the  next  four  sections,  automatic  procedures  for 
generating  variational  trial  electronic  wave  functions 
which  are  eigenf unctions  of  the  operators  ^^^^  (T^^ 
(if  J\_     =   0),  i,  ( if  the  molecule  is  homopolar),  ^   , 


and  $d        will  be  described. 
z 


4.   Construction  of  the  Slater  Determinants 

The  most  extensively-used  basis  functions  for 
constructing  the  MO ' s  (2.6)  are  the  atomic  orhiial  ( AO ^ 
functions  first  suggested  by  Slater  for  atomic  calcula- 
tions.  These  are  of  the  form 
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(4.2) 


(4.1)    (n,i,m)  =  (2n)    ^[(2n).']  ^  r'^-^e"^^  ^^m^^''^^ 

and  are  referred  to  as  Slater  orbltals  ( SO ' s )  or  Slater- 
type  orbitals  { STO ' s ) ,  depending  upon  whether  or  not  the 
[J, '  s  are  given  the  values  recommended  by  Slater  for  use  with 
atoms.   In  particular,  the  STO ' s  used  in  the  calculations 
reported  here  are,  in  spectroscopic  notation, 

,  3/2      -1/2  -|xr 

2s        =      [i^'^     {J>TT)      '       e 

o  —  5/2      -1/2  -\ir  a 

2p(r     =M.  TT     ^  e^cosS 

2p7r     =      [i-^^     (27r)      ^       e  ^      sm   9   e      ^    . 


These  five  functions  are  put  on  each  nucleus  with  the 
polar  axis  pointing  toward  the  other  nucleus.   This  gives 
a  total  of  ten  basis  functions  "X, , . . . ,  Xiq  which  are,  by 
the  manner  in  which  they  are  defined,  eigenfunctions  of  £    . 

C-i 

To  make  the  MO '  s  orthonormal,  the  overlap  matrix  _S  ,  whose 
elements  are 

(4.5)  S,j   =  /  XT  Ij  dx  , 

with   /  dT  denoting  integration  over  the  one-electron 
space  coordinates,  is  diagonalized  by  the  Jacobi 
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iteration-rotation  method  [20].   This  gives  an  orthonormal 
matrix  U   such  that 


4 

(4.4)  U   S  U   =   D 


is  diagonal.   Then  the  matrix 


(4.5)  T  -   U  D 


1 
2 


gives  a  transformation  to  a  set  of  MO ' s  expressed  as 
linear  combinations  of  atomic  orbltals  (LCAO), 


10 
(^.6)  ^.      -  TZ    ^i  t.  .  ,        j  =  1,2,. ..,10  , 


which  satisfy  the  orthonormality  condition 


(4.7)  J    i>l   4.J  dx   =   6^. 


If  the  molecule  is  homopolar,  the  orthonormalization 
procedure  can  be  made  to  yield  eigenfunctions  of  the 
operator  yC  .      Due  to  the  symmetry  of  the  AO's,  the  MO ' s 

so  constructed  are  also  eigenfunctions  of  Jl    .      Letting 

+ 

(T  and   tt   denote  the  MO '  s  with   1=0   and   ±1, 

z 

respectively,  the  MO '  s   (J)-,,  <t)p,  .  .  .  ,  ^^^   will  be  labelled 

(4.8)    IC-   10-  ,  2^,  2(r  ,  3(r,  J>cr  ,    tt",  tt",  tt^,    tt^ 
g    u'    g'    u'-^g'^u'   g'   u'  g'      u 
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In  the  usual  spectroscopic  notation.   When  a  calculation 
Is  completed,  the  numbers  1,  2,  and  3  on  the  resulting 
(J'-type  MO '  s  are  generally  assigned  In  the  order  of  the 
energies  associated  with  them,  starting  at  the  lowest.   To 
simplify  notation,  the  symbols  in  (4,8)  will  be  used  here 
whether  or  not  the  molecule  is  homopolar.   It  is  also  to 
be  understood  that  In  the  latter  case,  operations  having 
to  do  with  the  operator  ^  are  omitted. 

Multiplying  these  ten  MO ' s  by  a  and   P   spin 
functions  gives  a  total  of  twenty  MSO's,   ^-,  ,  ip^,...,    ip^Q, 
which  are  indexed  as  in  (2.7).   With  a  basis  for  a  one- 
electron  function  space  thus  established,  a  basis  for  an 

20' 


N-electron  function  space  is  taken  as  the  set  of   ,  ^ 

SD's, 

_1 
(^•9)      i^  =   (NJ)  ^  detj^^-^d),...,  ^^(N)]   , 


obtained  by  selecting  k.,  kp,...,k^   in  all  possible  ways 

from  among  the  integers   1,  2,...,  20.   To  avoid 

redundancy,  the  sets   ( k-,  ,  kp,...,k^)   will  be  ordered 

k-,  <  kp  <  ...  <  k^.   In  the  notation,   K  will  denote  such 

an  ordered  set  but,  when  convenient,   K  will  be  treated 

as  an  index   (  =  1,  2,.,.)   going  over  all  ordered  sets. 

As  a  result  of  the  orthonormallty  of  the  MSO's, 

Tp^,    ?//p,  .  .  .  ,    foQ'    ^^®   SD's    (4.9)    satisfy  the   orthonormallty 

condition 

-  31  - 


(4.10)        J     l^i^dV  =   6^L 


where   5„^   =   1   if  the  ordered  sets  K  and  L  are  the 
same  and   6^-r   =   0   if  they  are  not.   In  addition,  the 
SD '  s  are  already  eigenfunctions  of   >^  j   '^^^  ^^'■^    '^ 
whose  eigenvalues  can  easily  be  computed.   A  procedure 
for  forming  linear  combinations  of  the  SD ' s  which  are  also 
eigenfunctions  of  ^         and   CT   will  be  given  in  the 
next  three  sections. 

In  the  Hartree-Fock  SCF  method  for  closed  or  open 
shells,  one,  in  effect,  calculates  a  unitary  transformation 
R  of  the   T  matrix  (4.5), 


(4,11)  T^^^  =   T  R  , 


giving  a  set  of  MO ' s 

(4.12)  kf^   ^  YH    ^i  t?f 

SCF      ^  .       ,       OT.      tItSCF 


sue 


h  that  the  energy,   E    ,  for  a  single  SD  I^ 


constructed  of  the  MO ' s  (4.12),  attains  a  minimum.   Since 

,  SCF 
the  wave  function  IjJ     can  be  expressed  as  a  linear 

combination  of  the  basis  set  lj]_  ,   K  =  1,  2,...,  described 

K 

above,  the  wave  function  obtained  by  simply  finding  the 

best  linear  combination  of  the  l|j   's  will  be  better  than 

,SCF  ^ 

f     .   On  the  other  hand,  if  the  computational  effort 
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involved  with  such  a  large  number  of  SD ' s  is  too  formidable, 
or  if  SCF  MO ' s  are  already  available,  it  may  be  practical 
to  use  the  SCF  MO ' s  (4.12)  instead  of  the  MO • s  (4.6)  and 
restrict  the  set  of  SD ' s  to  those  which  contain  the  MO ' s 
having  the  lowest  eigenvalues  (orbital  energies)  of  the 
one-electron  operator  of  the  SCF  scheme.   In  this  case, 
the  automatic  procedures  described  here  for  constructing 
pure  state  functions  permit  the  use  of  far  more  SD ' s  than 
one  could  otherwise  include. 

In  the  following  methods  for  constructing  pure  state 
wave  functions,  the  only  properties  of  the  MO ' s  considered 
are,  of  course,  their  symmetry  in  the  molecule.   The 
author  would  like  to  suggest  that  the  use  of  some  other 
set  of  functions  having  the  correct  molecular  symmetry 
would  probably  be  better.   As  an  example,  MO ' s  of  the  form 

(4.13)         f(?,Ti)  exp  (-H?  +  ±m^) 
where 

i      =   (r^  +  r^)/R  ,     r]      =      ( r^  -  r^)/R 

not  only  seem  more  appropriate,  but  give  wave  functions 
quite  like  the  ones  used  by  James  and  Coolidge  [9]  and 
others  except  that  the  present  functions  would  not 
explicitly  contain  the  distance  between  pairs  of  electrons. 
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A  particularly  convenient  set  of  functions  of  the  form 
(4.13)  is  the  set  of  solutions  of   Hp  .   Such  a  choice  of 
basic  MO ' s  would  also  avoid  the  problem  of  "over- 
completeness"  which  results  from  the  fact  that  the  STO ' s 
form  a  complete  set  on  each  center.   The  effect  of  the 
overcompleteness  is  that,  with  a  large  number  of  STO ' s  on 
both  centers,  any  function  on  one  center  is  almost  a 
linear  combination  of  those  on  the  other  and,  as  a 
consequence,  there  is  a  large  indeterminacy  in  the  SD ' s 
which  becomes  more  serious  as   R,  the  internuclear 
distance,  decreases. 


5.   The  Projection  Operator  Technique 

The  use  of  projection  operators  by  the  methods 
developed  by  Lowdin  [21,  22]  not  only  is  conceptually 
simple  but  also  leads  to  systematic  procedures  which  are 
well-suited  for  automatic  machine  methods.   A  general 
description  of  the  fundamental  ideas  and  the  derivation 
of  some  necessary  formulas  are  given  in  this  section. 

Suppose  that  ^     is  an  arbitrary  function  and  that 
'^  is  a  Hermitian  operator  for  some  function  space 
containing  ^  .   Under  certain  conditions,  it  is  possible 
to  express  ^  in  the  form 

k   ^  ^ 

-  3J+  - 


where  _^  ,  ^  , . . „   are  elgenfunctions  of  '^  corresponding 
to  distinct  eigenvalues   A-,  ,  Ap,  o  .  .  .   When  this  Is  so, 
the  projection  operator 


(5.2)      Cr^      =       \     \     (A^,  -^)/(A^,  -A^)  , 
when  applied  to  ^  ,  yields 


(5.3)  (?!,$  =  c,^  . 

Thus,  an  elgenfunctlon  of  fyi     Is  formed  whose  eigenvalue 
Is  the  pre-selected   A,  .   The  denominator  In  (5.2)  gives  0^ 
(we  drop  the  subscript   k  for  the  rest  of  the  discussion) 
the  property 

(5.4)  ly'^   =  cr  . 

Suppose,  for  the  moment,  that  one  has  a  set  of 

orthonormal  functions  ill  ,  tIJ  ,  .  ,  ,   such  that  for  every  l|j 

-1  -2  -K 

in  the  set,  "^iIj    is  a  linear  combination  of  functions 
in  the  set.   Then,  if  0^    is  as  defined  in  (5.2)  one  can 
write 


(5.5)       ©I,  '    Crl^    '    Ii:   1^  c^L 

for  some  set  of  coefficients  ^vj-      The  elements  of  the 
overlap  matrix  of  the  projections   ©^   are 
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@^    0,.   dY     ^     J    iO"^    )    iO'l    )    dV 


K    ^L 


III     ^^      dV      =      ^  C^  ,^       /      III      llJ        dV 


(5.6)  =      C^^    . 


The  quantum-mechanical  "turnover  rule"  is  valid  here  since 
(^   is  Hermitian.   Since  the   C„^  's  are  the  elements  of 
the  overlap  matrix  of  the   0^.  's,  the  Schmidt  orthogonal- 
ization  procedure  can  be  applied  to  the   @^  ' s  by  simply 
triangularizing  the  matrix  of   C„^  's  by  fundamental  row 
operations.   This  process  discards  the  dependent  vs) „    's 
and  yields  a  new  set  of   ^^  's  and   C    's  having  the 
property  (5-6) .   Hence  the  formula 


Li   j  J-j 

(5.7)    -7^     ^L 

for  the  normalization  factor  Is  valid  both  before  and 
after  orthogonalization.   This  is  not  only  useful  as  a 
check  during  the  calculation  but  also  enables  the  program 
to  test  for  a  zero   fe)-.   during  the  orthogonalization  by 
simply  testing  the  L-th  diagonal  element  of  the   Cj,^ 
matrix. 
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It  will  be  convenient,  at  this  point,  to  introduce 
vector  and  matrix  notation.   A  set  of  functions  will  be 
written  as  a  row  vector  and  the  coefficients  of  a  linear 
combination  of  functions  will  be  written  as  a  column 
vector.   For  some  matrices,  a  symbolic  subscript  will  be 
used  to  denote  which  basis  is  used  to  define  the  elements, 
Accordingly,  (5.5)  and  (5.6)  can  be  written 

0  =   (y  I  =   I  C 
(5.8) 

y^  ®  ^  0  dv  -  c . 

Here,  "  ^    "  denotes  the  conjugate  transpose  and   /  dV 
indicates  that  the  integration  is  over  the  entire  space 
on  which  the  functions  are  defined. 
In  this  notation. 


-@    J   -   /^   -        J  - 


i'^!>/-(^idV  =   /  t^  ^1  dV  C  , 


where  the  properties  that  (^    is  Hermitian,  commutes  with 
^    ,  and  satisfies  (5-4)  have  been  used.   Hence,  the 
equation 
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(5.9)  H^=  HjC 

gives  the  relation  between  the  two  matrix  representations 
of  ^   .   Equation  {5.9)  Is  valid  If  .V  Is  replaced  by  any 
operator  which  commutes  with  ^9^  . 

The  Schmidt  orthogonallzatlon  procedure,  followed  by 
division  by  the  normalization  constant  ^'^vv    >    yields  a 
transformation  to  an  orthonormal  set  of  Independent  pure 
state  functions 

(5.10)  J  =  ®P  =  iC  P  =   I^  A 

where 

A  =   C  P  . 

Since  linearly  dependent  columns  of  _C   are  discarded,   _C 

is  a  rectangular  n  X  r  matrix,  where   n   is  the  number 

of  l|j   's  and   r   is  the  rank  of   C  .   The  matrix   P  , 
-K  -  - 

which  gives  the  transformation  to  orthonormal  functions, 
is  an   r  X  r   triangular  matrix  and  A   is  a  matrix  of  the 
same  dimensions  as   C_  .   From  (5' 9)  and  (5. 10),  one  gets 

H     =p''"h    P   =   p'^HCP. 

-5      --©-      — i"- 

Hence,  the  matrix  of  ^  with  respect  to  the  pure  state 
functions  can  be  expressed  in  terms  of  the  matrix  of  j{ 
with  respect  to  the  basis  IJJ  by  the  relation 
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5.11)  H    =   P  H   A 

-?     --llT- 


/^ 


The  SD's  of  the  previous  section  and  the  operators 
and   CT   have  the  properties  assumed  here  for  the 


basis  functions  IJJ   and  the  operator  "^  o   The  elements 

of  the  matrix  E         can  be  expressed  In  terms  of  the 

1 
Integrals 


o 

(5-12) 


^^(1)  H^  ^_g(l)  dv^ 


^'1  -  P 
'k  k'\        *,,   ^  .,*    /o^  /^      12 


£    £ 


z'  ±  —   r 

^k^i)  K'^^^{-F^)  n^^^  ^i'^2)  dv^  dv^ 


where   /  dv.   denotes  Integration  over  both  space  and  spin 
coordinates  of  electron   1  ,   Letting  K  =   ( k-,  ,  kg^.-.jk^), 
the  diagonal  elements  of  H    may  be  expressed 

i 

N      r/k.\  N  /k.    k.\ 

(5.13)  H^,     =    J-         ,         .    YZ.      (k!kO 

1=1       '-  1  J  =  l  +  1  1  J    ^ 

If   the    sets     K     and     L     differ  by   just   one   orbital,    k      in 
K      and      i      in      L    , 

(5-1*)  "kL     -      '-!)'"     [(^^^?=        Cyi 

k^A 
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where  s     is  the  number  of  Interchanges  required  in  the 
set   L   to  bring   i   to  the  position  of   k   in  the  set  K, 
If  K  and  L  differ  in  two  orbitals   k,k'   in  K   and 
i,i'   In  L, 

KL  ^ 

(5.15)  H,,  ^  (-i)^*^'"!]^:) 

KL 
where  s, -'^o/is  the  number  of  interchanges  in  L  necessary 

to  bring  £      and   i'   to  the  positions  of   k  and   k'  , 

respectively,  in  K.   If  K   and  L  differ  in  more  than 

two  orbitals,  then 


(5- 16)  Hj^^  =   0 


6,      Spin  and  (T       Projection  Operators 

Now  it  remains  to  apply  the  projection  operators  for 

the  spin  and  (if   L  =0)   (T       operators  to  determine  the 
'^  z         V   ^ 

matrices  A   and   P^  in  (5-11)  •   Of  the  various  methods 
for  forming  pure  spin  functions,  the  one  which  seems  most 
suitable  for  automatic  machine  computation  is  that  given 
by  Lowdin   and  used  here. 

Let  the  occupation  numbers, 

n^,  n^, . . .     ;         n^  =  0,  1,  or  2  ;    YIZ   n^  =  N  , 


See  Page  21,  Ref.  [21], 
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of  a  Slater  determinant  ( SD )  give  the  number  of  times  the 
respective  MO '  s  ^-.  ,    <!)„,...   appear  in  it  and  let  f  l|J  I 
denote  the  list  of  all  SD ' s  having  a  given  set  of  occupa- 
tion numbers  and  a  given  eigenvalue   s   of  ^a    •      The 
notation  [  K}  will  denote  the  list  of  ordered  sets  K 
corresponding  to  the  SB's  f  l|j  ]•  .   The  number  of   a 
electrons,   p.,  and  the  number  of   P   electrons,  v,  in  the 

singly  occupied  MO  '  s  of  the  SD '  s  in  \%'\       satisfy 

K 

(6.1)  2s   =  p.  -  V 

z 

(6.2)  N^   =   ii  +  V 


where  N   is  the  number  of  singly-occupied  MO ' s .   The 

s 

procedure  given  here  will  project,  from  \%\    ,    all  the 
eigenfunctions  of  xj         having   s  =  s  .   The  multiplicity 
of  the  corresponding  eigenvalue   s(s+l)   is 

(6.3)  M  =   2s  +1  . 

The  list  of  SD '  s  fijj  ^  ,  which  can  be  described  by  listing 
all  ways  of  distributing  \\.   a-spins  and   v  P-spins  over 

the  N    singly  occupied  MO's,  will  be  assumed  to  be  in 

s 

alphabetical  order,  i.e.,  the  order  in  which  the  list 

I  ijl  1-  would  appear  in  a  Greek  dictionary  if  just  the   a's 

and   P's  for  each  SD  were  written.   For  example,  if 


\\    - 


N   =   5  and  a  doublet  state,  i.e.,   M  =   2,  Is  desired, 
s 


then   s  =  -i  ,   M.   =   3  ,  and   v   =   2.   The  SD '  s  f  l|j  \ 

for  this  case  are  listed  in  the  first  column  of  Table  1. 

It  is  already  known  that  with  this  ordering,  the  spin 
projections  of  the  first 

(6.4)  ^  (""^ 


LL  +  1   , 


SD '  s  in  )  ill  ^  give  the  maximum  number  of  linearly 
independent  projections  of  the  whole  set.   There  are, 
therefore,  five  independent  doublet  states  In  the  present 
example.   The  projection  operator  which  projects  the 
s-state  component  out  of  a  determinantal  wave  function 

•X- 

having   s   =   s   can  be  written 
^   z 

(6.5)     ^Ja^^lP^]  =  nrrlZ  (-1)^  f^y'  [«^"P  P^l«^  P'"P]  > 

where  the  bra.cketed  expression  in  the   p-th  term  is  the 
sum  of  all  SD ' s  obtained  by  exchanging  p   a-spins  for 
P-spins  in  the  SD 

(6.6)  a. . .a  p. . .P   . 


Equation   All,  of  Ref.  [21] 
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TABLE  1.   INDEPENDENT  SPIN  PROJECTIONS  FOR  FIVE 
UNPAIRED  ELECTRONS  IN  A  DOUBLET  STATE 


-K 

&   aaapp 
s 

s 

s 

(J-  apaap 
s 

s 

aaaPP 

1/2 

-1/6 

-1/6 

-1/6 

-1/6 

aaPaP 

-1/6 

1/2 

-1/6 

-1/6 

1/6 

aapPa 

-1/6 

-1/6 

1/2 

1/6 

-1/6 

aPaaP 

-1/6 

-1/6 

1/6 

1/2 

-1/6  1 
1 

aPaPa 

-1/6 

1/6 

-1/6 

-1/6 

1/2 

aPPaa 

1/6 

-1/6 

-1/6 

-1/6 

-1/6 

paaaP 

-1/6 

-1/6 

1/6 

-1/6 

1/6 

PaaPa 

-1/6 

1/6 

-1/6 

1/6 

-1/6  ! 

PaPaa 

1/6 

-1/6 

-1/6 

1/6 

1/6 

PPaaa 

1/6 

1/6 

1/6 

-1/6 

-1/6 

The  elements  of  the   C_  matrix  for  the  spin  projection 
operator,  which  appear  as  the  coefficients  in  (6.5)  niay 
be  expressed  in  the  form 

•  1 


(6.7) 


C 


M 


KL 


|i,+l 


(-1) 


P 


KL 


P 


KL 


where   p^^   is  the  number  of  spin  interchanges  required 

KL 

to  change  JJ  into  l|J   .   Application  of  this  procedure 
K         L 
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to  the  first  five  SD ' s  of  the  list  in  the  example  gives 
the  matrix  of  coefficients  in  Table  1.   The  table  is 
arranged  so  that  the  K-th  column,  K  =  1,..,,    5>  gives 
the  spin  projection  of  the   K-th  SD.   If   L^  0,   C 
need  only  be  orthonormalized  to  give  the  matrices  A 
and  2  in  (5.11). 

Now,  suppose  that  a  rectangular  _C  matrix  is  formed 
by  the  above  procedure  and  that  L   =0.   Then,  one  must 
apply  the  projection  operator 


(6.8)         ^^  =  |(i  ±  <y^\ 


+ 


where  +  is  used  for  a  ) state  and   -   is  used  for  a 

)      state.   It  is  desirable  to  express  the  effect  of 

(y       in  terms  of  the  elements  of   C_  ,  i.e.,  without 
requiring  coefficients  of  the  linearly  dependent  spin 
projections.   The  effect  of   (T   on  the  SD '  s  I  IJJ  ]   is 
to  replace  the  occupied  MO ' s  by  their  complex  conjugates. 
Due  to  the  convention,  adopted  here,  of  keeping  the  indices 
of  both  the  MO ' s  and  the  MSO ' s  in  numerical  order  (2.7), 
the  sign  change,  which  may  result  from  putting  the 
conjugates  of  the  MSO ' s  in  this  order  after  operating  with 

(T  ,  is  the  same  for  the  whole  set  {ijj  |  .   The  effect 


of   CT   can  be  expressed 
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(6-9)      ^v{y  =  {i;i  ^^[SkJ- 

This  relation  defines  a  one-one  correspondence  between 
the  list  of  ordered  sets  i  K  j  and  the  list  of  ordered 
sets  [K'I  ,   Let   [cT  |   be  defined  by 

(6.10)        la-J  [t]     -    [t  ,]  . 


■K-'     ^  -K  ' 


2 
Since  rr   =   1, 


(6.11)       Icr^l  fi|J  I  =  ft 


Therefore,  when   (5^   operates  on  the  SD '  s  In  [  ]^  I  ,  it 
can  be  expressed 

(6.12)        Cr^     =  1(1  +  sgn  Icr^l  ) 

where   sgn  =   ±1   is  the  product  of  the  alternative 
signs  in  (6,8)  and  (6,9) •   Now  let 


(6.13)         ^=   ^s  ^v  =   ^v  ^s  • 

The  second  equality  holds  because  X)    and   (T   commute 

If  the  operation   (T"   changes  the  occupation 
numbers  of  j''^    X  ,    then  the  lists  (  ]j[  }  and  f  l|_  \     have 
no  SD's  in  common.   The  result  of  applying  (y       to  the 
spin  function 
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(6.14)         ^  i^  =  _YZ^  i^  c^j^ 


^  ~K     l1^^  -^ 


IS 


6.15)    (Ti,   =   ZZ   (1  I  ^LK  +  t  ,  SS^  -I  ^LK^ 


On  the  other  hand,  if   0"  does  not  change  the 
occupation  numbers,  the  lists  /  1^  }   and  j"  l|j  1   have 

the  same  SB's  in,  perhaps,  a  different  order.   In  this 
case,  (6.14)  can  also  be  written. 


^      L'£-^K'J-   ^ 

Using  both    (6.14)    and    (6.l6),    and   the   relation    (6.11), 
one   gets 

/5't|j      =    ^  f  }        I    c.^  +  sgn  Icr  I     J  ^      a  ,^ 

^-K  2|l%3-L      LK  ^      I     v'    ^4^,^    ^,      L'K 

=      If   IZ     1    C        +    sgn     YIZ       1    C 
Therefore, 

(6.17)      (^1        =     m     i     |(C^j^   +    sgn   C^.j^)       . 
K  Lf-TLj      L 

The  coefficients  of  the  SD ' s  in  (6.I5)  or  (6. 17),  as 
the  case  may  be,  give  the  C     matrix  of  (5-9)  for  the 
operator  (y  .      When  calculating  the  matrices   A   and   P 
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of  (5- 11)  by  the  orthogonallzation  procedure,  the  dependent 
projections  resulting  from  the  application  of  Qf"       are 
found  and  discarded. 


7.   The  Machine  Program  For  Constructing  Pure  State  Wave 
Functions 
The  I.B.M.  704,  for  which  these  calculations  have 
been  programmed,  has  certain  Boolean  operations  which  treat 
the  36-blt  binary  words  as  sets.   That  Is,  a   1   or  0   In 
bit  position  n   can  mean  that  object  number  n   Is  or  Is 
not  In  the  set.   Letting   A  and  B  be  two  sets  represen- 
ted In  this  fashion,  a  single  machine  operation  can  compute 
their  Intersection, 

C   -   A  *  B  , 

which  Is  a  word  having   I's  where  both  A   and  B  have 
I's  and   O's  elsewhere.   Another  operation  gives  the  union 
of   A   and  B, 

C   =   A  +  B  , 

which  Is  a  word  having   I's  where  either  A   or  B  have 
i's  and  O's  elsewhere.   The  complement  of   A, 

C  =   (-A)  , 
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Is  formed  by  changing   O's  to   I's,  and  vice-versa,  in  A. 
The  notation  used  here  for  the  Boolean  operations  is  that 
used  in  programming  them  in  the  latest  FORTRAN  coding 
system. 

Since  the  Slater  determinants  described  above  are 
ordered  sets  of  objects,  no  one  of  which  can  appear  more 
than  once,  it  is  convenient  to  represent  them  by  binary 
words  with   0   or   1   representing  the  absence  of  presence 
of  a  function  in  the  determinant.   The  Boolean  operations 
then  enable  one  to  perform  tests  and  operations  on  these 
sets.   Specifically,  an  SD  is  represented  in  the  first 
twenty  bits  of  a  word  with  bit  positions  corresponding  to 
the  twenty  one-electron  functions. 


l<ra,  lo^p,  icr^a,  l(/^p,  2(ra,  2(rp,  2cr^a, 

(7.1)      20;p,  3Cga,  J>(r^^,  J>6-^^>  3Cr^^,  TT^a,   tt^P, 

—      —      +  +     +  + 

U  '     U  '     g  '  g  '     U  '  U  ' 


in  that  order.   For  example,  the  bit  pattern 

A  =   1011010  ...  0 

represents  the  determinant 

_1 

A     =      ik:)    ^  det    fKTa,    ICTa,    ICT  p,    2Cr  p  ]• 

L      g    '         u    '         u    '         g     -^ 
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To  illustrate  how  tests  and  operations  are  performed, 
consider  the  formation  of  |(r  |A,  the  complex  conjugate, 

without  sign,  of  a  determinant  A.   The  operation  JC  | 

±      + 
simply  replaces  each  v        by  tt   and  can  be  performed  by 

shifting  the  contents  of  the  four  w        bit  positions  to 

the  four  tt   positions,  and  vice-versao   To  observe  the 

the  specific  steps  involved,  let   TS,  TM,  and  TP  be 

words  having   I's  in  the  (T,  tt",  and  tt  positions, 

respectively,  and   O's  elsewhere.   Then,  using  Boolean 

operations,  the  machine  forms 

AS  =   A  *  TS 
(7.2)  AM  =   A  *  TP   shifted  left  4 

AP  =   A  *  TM   shifted  right  4   . 


Then   C  A   is 
'  v' 


(TO)  AS  +  AM  +  AP  . 

The  sign  in  equation  (6.9)  is  Just 

(^^4)  (_^)(92(AM)^(AP)) 

where  '[    is  a  function  which  computes  the  number  of   I's 
in  a  word.   In  the  same  notation,  the  calculation  of  the 
eigenvalue  of  cX.         can  be  described  by  the  formula 
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(7.5)  L^   =  '^(AP)  -  ^(AM)  . 

The  program  for  computing  the  pure  state  wave 
functions  accepts  as  Input,  N,  the  number  of  electrons, 
M,  the  multiplicity,  and  the  eigenvalues  of  the  three 
operators  <!,      ,    -^-'  ,  and   u""'  .   The  multiplicity  determines 
the  eigenvalues  of  XJ    and  ^  .   The  number  of 


occupied  orbltals,  N  ,  can  assume  all  Integral  values  such 
that 


(7-6)  I  ^  N  <_   Mln  (10, N) 

where   I   is  the  Integral  part  of  -^(N+1).   For  each  N  , 
which  can  be  considered  fixed  for  the  rest  of  the  present 
discussion,  the  program  calculates  the  number  of  doubly 
and  singly  occupied  orbltals. 


and 


d         o 


N    =   N   -  N ,  , 
s      o    d  ' 


respectively.   The  numbers  of   a  and   P   electrons  in 
the  singly  occupied  MO ' s  are 


and 


M.   -   ^(Ng+M-l 


s   ^  ' 
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respectively.   To  form  the  bit  patterns  for  the  SD's,  a 
subroutine  was  written  which  accepts  a  given  M-,   and  M^ 

and  generates  a  list  of  (  ^^  )  binary  words;  one  word 

for  each  way  of  distributing  M-,   I's  among  the  first  M^ 
bit  positions.   The  list  would  be  in  alphabetical  order 
if   1   and   0   were   a  and   P,  respectively.   This  sub- 
routine is  used  to  generate  three  lists,  D,  S,  and  A 
with  elements 


(7.7) 


^i ' 

1   =    1,    2,...,   (^J 

J    =    1,    2,...,^^ 

s 

\  ' 

/     N         N 

k  =    1,    2,...,l       ^ 

The  list   D  gives  all  ways  of  selecting  N,   doubly 
occupied  MO ' s  from  the  ten  available  ones,  S  gives  all 

ways  of  selecting  N    singly  occupied  MO ' s  from  among 

s 

the   10-N^   remaining  MO's,  and  A  gives  all  ways  of 

distributing  [i   a-spins  over  the  N    singly  occupied 

s 

MO's.   The  triplets  D. ,  S.,  A  ,  with   i,  j,  k,  going 
over  the  ranges  in  (7.7)  describe  all  SD ' s  having  N 
occupied  MO's. 
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The  _C   matrix  for  the  spin  projection  operator, 

which  will  be  denoted  by   C  ,  depends  only  on  A;  not  on 

s 

D  and   S.   The  elements  of   C   are  calculated  by  using 

s 

(6.7)  where   p,^^   Is  simply  the  number  of   1 '  s  In  the 
result  of  the  Boolean  operation 

on  the  spin  patterns   A,   and   A.   of  K  and  L.   This 
Is  simply  a  way  of  counting  the  number  of  positions  where 
A,   has  a     and  A.   has   p.   If  L   7^  0,  then  C    Is 

K  i/  Z  — S 

orthonormallzed,  yielding  the  matrices   A  and  _P  of  (5. 11) 
for  an,  as  yet,  unspecified  set  of  MO's.   Otherwise,  _C 
Is  saved. 

Then,  for  each   i,J   In  the  ranges  defined  by  (7-7), 
D.   and   S.   are  examined  to  see  If  they  give  SB's  having 


the  required  eigenvalues  of  ^    and  ^  .   If  not,  the 


of^ 


/N, 


combines  D. ,  S.   with  A   for  k  =  1,  2, . . . ,       ,  to 


program  goes  on  to  the  next   l,j.   Otherwise,  a  subroutine 

fs^ 

In- 
form a  list  -T  ^  I  in  terms  of  bit  patterns  according  to 

the  scheme  in  (7.I).   If  L   7^  0,  the   A  and   P 

matrices,  which  were  calculated  for  the   C   matrix  and 

— s 

saved,  give  the  result  in  the  form  required  by  (5>11). 
Otherwise,  the  (y  operator  is  applied  by  the  use  of 
equations  (6.I5)  or  (6.17).   This  Involves  simple  row 
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operations  on  _C   and  the  formation  of  the  complex 
conjugates  y  ill   7  by  the  operations  described  in  (7-2) 
and  (7' 3)'   The   C_  matrix  resulting  from  these  operations 
on  C_        Is  then  orthonormallzed,  dependent  columns  are 
discarded,  and  the   A  and   P_  matrices  required  by  (5.II) 
are  calculated. 

The  list  ■[  i  ]  ,  together  with  the  list  j   f  ^']   if 
L   =   0,  will  be  called  a  "configuration".   From  the 
above,  it  is  seen  that  the  result  of  operating  on  any  SD 
of  a  configuration  with  (^  is  a  linear  combination  of  SD ' s 
in  that  configuration  only.   Therefore,  the   A   and   P^ 
matrices  of  (5-11)  can  be  partitioned  by  the  configurations 
so  that  they  appear  as  block  diagonal  matrices.   Treatment 
of  them  as  such  produces  a  great  saving  in  storage  and 
calculation. 

The  elements  of  H.!.   in  (5,11)  are  calculated  by 
subroutines  which  use  the  Boolean  operations  to  examine 
the  bit  patterns  for  the  l|j   '  s  in  the  book-keeping 
process  which  selects  the  one-  and  two-electron  integrals 
needed  in  equations  (5"13)j  (5'1^);  and  (5' 15)" 

In  order  to  perform  the  elimination  of  dependent 
projections,  it  is  necessary  to  detect  the  presence  of 
identically-zero  matrix  elements  during  the  calculation. 
This  means  that  an  exact  representation  of  all  numbers 
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must  be  carried  throughout  the  calculation.   To  do  this, 
the  matrix  elements,  which  are  rational  numbers   a/b   are 
represented  In  the  machine  as  pairs  of  Integers   (a,b) 
and  the  calculations  are  performed  In  Integer  arithmetic  by 
FORTRAN  subroutines  written  for  the  purpose. 

8.   Averaging  Over  Electronic  Coordinates 

Nearly  all  important  properties  and  perturbation 
theory  estimates  of  weak  Interactions  are  given  at  each 
fixed  R   in  terms  of  integrals  of  the  form 

(8.1)  P(R)   =   r  I*  (P  t'  dV 

where  (p  is  a  sum  of  mono-electronic  operators 

N 

(8.2)  (P  =   IZ  (P(i) 

1=1 

and  JJ   and  1^   are  two  wave  functions  for  the  system. 

If  IJJ  =  £  ,  then  (8.1)  gives  the  average  value  of 
^O    when  the  system  Is  in  a  situation  characterized  by  the 
normalized  wave  function  f  .   If  ^     ^     IJJ  ,  (8.1)  gives 
Information  about  the  system  in  a  transition  from  the 
state  %_     to  the  state  %_   .      The  integral  (8.1)  for  any 
operator  of  the  form  (8.2)  can  be  obtained  in  terms  of 
the  function 
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(8.3)   p(l',l)  =    n    f    JJ*(1',2,.  .  .,N)l|j'(l,2,.  ..,N)  dV^^^  , 

where   /  dV     denotes  integration  over  all  electronic 
coordinates  except  those  of  electron   1.   If  1^  =  ^  , 
(8.3)  is  called  the  "density  matrix"  and  if  I^  ^     t  ,  it 
is  called  the  "transition  matrix".   In  this  section,  we 
need  not  distinguish  between  the  two.   If   1'   =   1,  and 
jj  =   IJJ  ,  (8.3)  gives  the  density  of  electronic  charge  at 
the  point  corresponding  to  the  space-spin  coordinates   1. 
The  primed  one-electron  coordinates  are  used  to  allow  for 
operators  on  the  phase  space  of  the  system. 

The  integral  (8.1)  can  then  be  given  in  terms  of 
p(l',l)   by 


(8.4)  /  (p  (1)  p(l',l)  dv 


1  ' 


where  the  notation  implies  that  (P  ( 1 )   is  to  operate  on 
only  the  unprimed  coordinates  of   p(l',l)   after  which  the 
primes  are  removed  and  the  integration  is  performed. 

To  get   p(l',l)   in  terms  of  the  MSO '  s  i/-^,    f^,..., 
consider  the  form  of  the  wave  function  involved. 


(8.5)  i  -  y~  c  ?   . 

K       K 

The  ^   's  are  the  pure  state  functions  in  (5. 10), 
-K 
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(8.6) 


so  that  (8.5)  can  be  written 


(8.7) 
where 


t 


L^L 


(8.8) 


B. 


K 


^LK  ^K 


Hence,    (8.3)    can  be   written 


9)        pdM 


^5^kB. 


-K  -L 


In  order  to  evaluate  (8.9),  one  needs  the  following 
relation: 


det  Iflil),.  ..  ,ip^{m)l    det  |t//-[(1),.  .  .,^j^(M)]  dv^...d^ 


M 


=  Mi  /  ^*(1),...,  ^*(M)det[^|(l),...,  ilf^{n)'}dv^...dv^ 


=  MJ  /  det  [V/i(l)^i(l),..-,%(M)^j^(M)(  dv^...dvj^ 


(8.10)   =  M. 


^11  • • •  ^IM 


^Ml  . . .  ^MM 


where 
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(8.11) 


k^ 


^,  (1)  ^,(1)  dv 


k 


I, 


1 


Expanding  the  determinants  In  (8.9)  in  terms  of  their 
first  rows,  one  obtains, 


p(l',l)  =  [(N-1).']-^  YH  \   ^L  ^.        ^k  ^^'  ^  ^^  ^^^ 

K,L       i,J=l    i 


(8.12) 


X  J   d^(k.)  d^(..)  dv^i) 


where  K  =  (k-,,...,k^)   and  L  =  (i,,...,i^)   and  where 
dT^(k.  )   is   (-1)      times  the  first  minor  of  the  element 


^^  (1)   in  the  SD  I^  .   By  the  relation  (8,10),  it  is  seen 

i 
that 


~K 


(8.13) 


d^(k.)  dL(ij)  dv^.-.dv^ 


1  J 


where  D 


KL 

k.^  . 
1  J 


is   (-1 


li+J 


times  the  first  minor  of  the 


k.,^.   element  of  the  determinant  of  the  "overlap"  matrix 
1   J 


(8.1if) 


S   = 


\h 


Vi 


.     s 


^l^N 


^N-^N 


whose  elements  are  defined  by  (8,11).  Hence, 


N 


(8.15)   pdM)  =  ]^W    H 

K,L       1,J  =  1   "l"j   '"i 


^k^.  <.'!■'  "^.'l 
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Changing  the  order  of  summation  in  (8.I5)  gives 

(8.16)        pdM)    =   y~  r^^j,  i^lil')  f^n) 

k,  £ 
where 

I 1       ITT 

^^     KeTk)  LTU)      ^  ^   ^^ 

The  summation  is  over  all  K   containing   k  and  all   L 
containing  £. 

The  integral  (8.1)  can  then  be  expressed  in  the  form 

(8.18)     p(R)  =  ^  rki/<(i')(r^(i)  ^>)  d^i 

in  terms  of  one-electron  integrals.   The  sum  in  (8. 17)  can 

be  formed  by  going  over  all  pairs  of  determinants  t      ,    f_ 

K     L 

and  including  the  K,L   term  of  (8. 17)  in  the  sum  for  each 

P,  „   for  which  ke  K   and   ieL  .   In  the  scheme  for 
'  kX 

representing  SB's  as  bit  patterns,  as  described  in  Section 
7,  a  simple  Boolean   "  *   "   operation  is  used  to  test  for 
the  presence  of  an  MO  ^,   in  a  determinant  1^   . 

If  the  one-electron  MO '  s  used  in  constructing  IJJ 
belong  to  the  same  orthonormal  set  as  those  used  in 
constructing  ijl  ,  then 


(8.19)  \^     =     \, 


and,  therefore. 
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KL 


where  K-k  and  L-£      denote  the  sets  K  with  k  omitted 
and  L  with  £      omitted  and  s,  .  is  defined  In  (5.14). 

In  the  special  case  of  properties  for  which  (y 
commutes  with  the  operators  whose  projection  operators 
were  used  in  obtaining  the  ^   's,  a  simplification  is 
available  which  amounts  to  allowing  one  to  use   C 


K 
Instead  of  B^   in  (8.17)  to  calculate  a  set  of  l-i^p 


's 


which  can  be  used  in  (8.l8)  to  calculate  the  average  value 
of  (p. 

9.   Optimization  of  Non-Linear  Parameters 

The  calculation  of  optimal  linear  parameters  of  a 
trial  wave  function  for  a  given  pure  state  consists  in 
solving  the  matrix  eigenvalue  equation 

(9.1)        LI  ^\l  -   \l  e)  ^l  =  0 


where  H^-j-   is  the  matrix  of  the  Hamiltonlan  of  the 
system  with  respect  to  a  set  of  orthonormal  pure  state 
functions  ^^  ,  K  =  1,  2,.,.  .   The  eigenvalues, 
E„  ^  E-,  ^  Ep  j^  ,  .  .  ,  and  their  corresponding  eigen- 
functions. 
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9.2)         i!        =  YZ  ^j"^   ^r  '  V  =  1,  2, 


are  approximations  to  the  eigenvalues  and  elgenfunctions 

of  the  wave  equation  (l.-l)  for  the  pure  state  under 

consideration.   Furthermore,  according  to  the  general 

separation  theorem   [25],   E^   is  an  upper  bound  of  the 

eigenvalue  it  approximates.   This  procedure  optimizes  the 

linear  parameters   C-.   of  the  wave  function,  but  one  must 

also  consider  non-linear  parameters  such  as  the  orbital 

exponents  of  the  AO's. 

Letting  [i      represent  any  set  of  non-linear  parameters 

in  the  wave  function,  the  eigenvalues  of  (9-1)  may  be 

regarded  as  functions   E  (|j.)   of  \i .      Then,  since   E  (ij.) 

is  an  upper  bound  to  the   v-th  exact  eigenvalue,  one  can 

improve   E  (|-l),  for  some  particular  value  of   v ,  by  finding 

the  point  (i   =   p.   where  E  (^i)   attains  a  minimum. 

The  corresponding  wave  functions,   l|j  (|j,),   v  =  1,  2,.,., 

~v 

— 

A  rough  outline  of  a  proof  of  the  separation  theorem  can 
be  formulated  by  regarding  Equation  (4,1)  of  Chapter  IV, 
below,  as  the  matrix  formulation  of  (9.1).   The  results 
In  Chapter  IV  show  that  if  the  basis  functions  for  (9-1) 
are  taken  from  a  complete  set  of  functions,  the  addition 
of  one  more  function  gives  a  new  set  of  eigenvalues  such 
that,  except  for  the  maximum  one,  each  is  lower  than  or 
equal  to  an  eigenvalue  of  the  previous  set.   Therefore, 
the  addition  of  functions  of  the  complete  set,  one  by 
one,  yields  a  non-increasing  sequence  of  approximations 
to  each  exact  eigenvalue. 
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are  orthogonal  to  each  other  if  all  have  the  same  [i ,    but 
the  optimized  wave  functions  IJJ  ill    )      will,  in  general, 
not  have  this  property.   This  means  that  in  the  calculation 
of  transition  matrix  elements  and  perturbations  one  will, 
at  least  in  effect,  use  wave  functions  which  are  orthogonal, 
but  less  accurate  than  the  l|j  (|j,  )  's.   In  the  case  of 
perturbation  calculations,  the  favorable  aspect  of  such  a 
procedure  is  that  one  can  select  and  optimize  a  particular 
eigenstate  for  the  unperturbed  system  and,  if  the  perturba- 
tion is  small,  the  effect  of  the  non-orthogonality  will 
also  be  small.   However,  the  projection  operators  do  not 
depend  upon  \i ,    so  the  wave  functions  IjJ  (|-t  )   will  be 
orthogonal  to  the  eigenfunctions  for  other  pure  states. 

The  calculation  of  electronic  wave  functions  for  even 
a  single  set  of  non-linear  parameters  requires  a  consider- 
able amount  of  calculation  and  one  will,  in  practical 
applications,  be  quite  limited  in  the  number  of  such 
parameters  one  can  vary  and  in  the  number  of  eigenvalues 
one  can  optimize.   Thus,  the  convenience  of  using  the 
small  number  of  basis  functions  in  the  calculations 
described  here  is  somewhat  lessened  by  the  fact  that  the 
orbital  exponents  are  quite  critical  in  the  determination 
of  energies  and  other  molecular  properties.   This 
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necessitated  putting  a  somewhat  large  amount  of  computa- 
tional effort  Into  the  procedure,  described  In  the  next 
section,  for  minimizing  a  function  of  several  variables. 

10.   Numerical  Method  for  Minimizing  a  Function  of   n 

Variables 

In  numerical  work,  one  frequently  encounters  the 
problem  of  having  to  determine  a  point  where  a  given 
function  F(v),   v  =  ( v-,  ,  Vp,  ,  ,  .  ,  v  ),  whose  values  can 
be  calculated,  attains  a  minimum.   The  amount  of  computing 
effort  Involved  in  obtaining  a  detailed  knowledge  of  the 
behavior  of  a  function  of  several  variables  from  its 
numerical  values  can  be  formidable  if  one  does  not  employ 
an  effective  means  of  utilizing  the  information  obtained 
in  successive  evaluations  of  the  function.   The  computing 
machine  subroutine  described  in  this  section  uses  a 
process  of  fairly  general  applicability  for  finding  a 
relative  minimum  point  of  a  given  function.   The  procedure 
is  best  described  by  first  giving  a  brief  outline  of  its 
essential  features.   On  each  iteration,  a  quadratic 
function  ,   G(v),  coincident  with  F(v)   at  the  best 


Some  basic  ideas  for  fitting  and  using  quadratic  surfaces 
for  minimizing  a  function  of  several  variables  are 
discussed  in  Ref,  [2^]. 
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(10.1)  m  =  -|(n+l)(n+2) 

calculated  values  of  F(v),  is  obtained  by  determining 
the   m  coefficients  of 

(10.2)  G(v)   =   a  +  b'^v+-|v'^Cv, 

where   b_  and   v;  are  column  vectors,  C   is  a  symmetric 

matrix,  and   "  +  "   denotes  the  transpose.   The  symmetric 

matrix  C_     is  dlagonalized,  giving  the  eigenvalues   C.  , 

i  =  1,  2,  .  .  .  ,  n  of   C_  .   If  the  points  used  to  determine 

(10.2)  are  in  a  sufficiently  small  neighborhood  of  a 

minimum  point  of  F(v),  then   C.  >  0   for   i  =  1,  2,...,  n 

and  the  minimum  point  of  (10.2)  gives  a  good  approximation 

to  the  minimum  point  of   F(_v)  .   If  the  condition   C.  ^  0, 

i  =  1,  2,..,,  n  does  not  hold,  a  new  point  is  found  by 

going  from  the  best  calculated  point  in  the  direction  of 

the  negative  gradient  of   G(v)   at  that  point. 

To  start  the  iterative  procedure,  the  program  is 

given  an  initial  point  _v   and  positive  step-sizes, 

d-,,...,d  ,  for  all  variables.   Then  F(v)   is  calculated 

at   v^  =  ( v-i  ,...,v   )   and  at  the   2n  points 
— o     lo     '  no  ^ 

(v^^,...,  v^^  ±  d^,...,  v^^),  i  =  1,  2,...,  n.   From  these 
computed  values. 


-  63 


S.  =  -  sign  [F(v,  ,...,v.   +  d.,...,  v   ) 
1       ^      lo'    '  lo    i'    'no 

-  F(  V-,  ,  ....  V.   -  d.  ,  .  .  .  ,  V   )  1 
lo'    '10    i'    '   no  ■' 

is  determined  and  F(_v)   is  calculated  at  the  -pnln-l) 

points  (v.,  ,...,v.   +  S.d.,...,v.   +  S.d.,...,v   ).   This 
^        lo'    '  lo    1  i'      JO    J  J      no 

provides  the   m  points  needed  to  fit  the  first  quadratic 

surface. 

On  each  iteration,  the  program  has  an  ordered  list 

of  points   V.   and  calculated  values   G.   =   F(v.)  , 
-J  J       -J 

j  =  1,  2,...,    m',  where   m'  =  m  or  m+1,  such  that 

G-,  <  G„  *^  ...  <  G  ,  .   The   m  independent  coefficients 
1  —  2  —     —  m'  ^ 

of  (10.2)  are  determined  by  the   m'   equations 


(10.5)  G   =  G(v  )  ,  j  =  1,  2,...,  m' 

J       J 


As  it  will  be  shown,  the  rank  of  the  system  of  equations 
(10,3)   is  at  least   m  at  the  start  of  each  iteration. 

Equations  (10. 3)  are  solved  by  the  Gauss  elimination 
method,  taking  them  in  the  order   J  =  1,  2,...  .   By  this 
method,  the  vanishing  of  the  coefficients  of  the   i-th 
equation  during  the  elimination  process  implies  that  the 
i-th  equation  is  redundant  or  inconsistent.   In  this  case, 
the   i-th  point  is  discarded  from  the  list  and  the 
H+1,    i+2, . . . ,  m'   points  in  the  list  and  the  corresponding 
equations  in  (10. 3)  are  re-numbered   i,  i+1,...,  m'-l.   If 
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m'  =  m+1   and  the  first   m  equations  of  (10. 5)  are  of 
rank  m,  then  the   (m+l)-st  equation  is  simply  not  used. 
A  necessary  and  sufficient  condition  that   G(_v)   have  a 
minimum  at  a  point  is  that  the  gradient 

(10.4)         V  G(v)   =   b  +  C  V 

be  zero  and  that  the  eigenvalues  of   C   satisfy 


(10.5)  c'  >  0  ,  i  =  1,  2,  .  .  .  ,  n  , 


at  the  point.   Therefore,  the  matrix  C      is  diagonalized, 
yielding  an  orthonormal  matrix  R,  whose  columns  are  the 
eigenvectors  of  C,    and  a  diagonal  matrix   C_'   whose   i-th 
diagonal  element  is  the  eigenvalue  of   C_  corresponding 
to  the   i-th  column  of   R.   Hence 

(10.6)  C   =   R  C'  R^  . 

If  all  of  the  eigenvalues  of  _C   are  positive,  then 
the  minimum  point, 

(10.7)  V  =  -  C"-^  b  =  -  R  C'"^  R^  b   , 

of  G(_v)  is  calculated.  If  _v  is  within  a  distance  d, 
to  be  defined  later,  of  the  best  calculated  point,  _v-,  , 
then   v'  ,  the  next  point  at  which  to  calculate  F(v),  is 
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set  equal  to   v^  .   Otherwise,  the  gradient  (10.4)  Is 
evaluated  at  _v,   and  _v'   Is  taken  as  the  point  at  a 
distance  d  from  v_,   down  the  gradient  VG{_v-,).   Thus, 

(10.8)      V'   =   v^  -  dVG(v-j_)/|  VG(v-j_)  I  , 

where   |VG(v_^)|   denotes  the  length  of  the  gradient. 

On  the  first  Iteration,  and  on  each  Iteration  for 
which  the  previous  one  gave  an  Improvement,  the  parameter 
d   Is  the  distance  between  the  two  best  calculated  points. 
If  the  previous  Iteration  did  not  yield  an  Improvement 
and  If  _v'   of  the  previous  Iteration  was  obtained  by 
going  down  a  gradient.  It  Is  evident  that   d  was  too 
large  and  the  new  d   Is  taken  as  half  the   d   of  the 
previous  Iteration.   In  the  only  remaining  situation, 
where   v'   of  the  previous  Iteration  was  obtained  as  the 
minimum  point  of  a  quadratic  surface  and  did  not  yield  an 
improvement,  the  old  d   is  used. 

Then,   F(_v'  )   is  calculated  and  the  test  for 

convergence,  as  described  below,  is  applied.   If  the 

convergence  criteria  are  not  satisfied  and  if  F(  v_'  )  ^   G  , 

the  new  point   v_'   and  value  F(_v'  )   are  placed  in  the 

appropriate  place  in  the  list  of  calculated  values.   All 

values  G.   in  the  list  for  which  F(_v' )  <  G.  are  moved 
J  J 

upward  one  position  and  renumbered.   As  a  result  of  the 
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method  of  solving  (10. 3),  the  m  equations  (10. 3),  before 

the  insertion,  are  of  rank  m.   Therefore,  the  rank  after 

Insertion  is  at  least  m.   On  the  other  hand,  if 

F(v')  >  G  ,  the  new  value  is  not  put  in  the  list.   Instead, 
—      m 

the  calculation  of  a  new  quadratic  surface  is  by-passed, 
d   is  replaced  by  -pd  ,  and  the  program  goes  back  to  the 
point  where  it  uses  the  gradient  to  find  a  new  point.   In 
other  words,  it  goes  down  the  same  gradient  as  on  the 
previous  step,  but  only  half  as  far. 

The  usual  experience  is  that  one  starts  with  points 
far  apart  to  get  the  maximum  searching  effect  of  the 
program.   Therefore,  quadratic  surfaces  fitted  in  the 
first  iterations  give  a  poor  description  of  the  behavior 
of  F(v).   Then,  as  the  computed  points  begin  to  cluster 
about  the  region  of  a  minimum,  a  fairly  good  description 
of  F(_v)   about  the  minimum  is  given  by  the  quadratic 
surface,  convergence  becomes  of  second  order,  and  the 
points  rapidly  come  so  close  together  that  significance 
in  the  calculated  coefficients  of  the  quadratic  surface  is 
lost.   Hence,  the  quadratic  surface  itself  is  not  very 
useful  in  forming  convergence  criteria. 

Among  the  arguments  of  the  computer  subroutine,  are 
two  parameters,   NP   and   e  .   If  NP  >  0,  the  subroutine 
will  exit  when  the  best  NP   calculated  values  of  F(v) 
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lie  within  e      of  each  other.   If  NP  =  0,  the  subroutine 
consults  a  sense  switch  to  see  if  it  should  exit.   This 
allows  the  user  to  use  his  own  judgement,  during  a 
calculation,  to  decide  when  the  procedure  has  converged. 
Then,  if  further  calculations  are  performed  with  similar 
functions,  the  experience  obtained  should  enable  one  to 
select  an  NP   and  an   e  . 
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Ill   THE  CALCULATION  OF  MOLECULAR  INTEGRALS 

1.   Integrals  of  the  STO ' s 

By  performing  the  trivial  integration  over  spin 
coordinates,  one  can  obtain  all  of  the  integrals  involved 
in  the  calculation  of  electronic  wave  functions, 
molecular  properties,  and  perturbations  in  terms  of 
integrals  of  the  MO ' s 

(1.1)  i>,    =  5Z  TCfcC,, 

k 

which  can,  in  turn,  be  expressed  in  terms  of  Integrals  of 
the  basic  set  of  STO ' s  (II  4,1)  on  the  two  centers.   For 
a  one-electron  operator  (P  ( 1 ) , 

(1.2)    r  <t)!(l)(P(l)  (l).(l)  dT,   =  5    C  .  C  .  P 

J   ^1   ^      ^J       1     ^-—   pi   qj   pq 

where 

(1.5)         Ppq  =  /7;(i)(p(i)  Xqd)  dT^ . 

Here,  dx-,   denotes  integration  over  the  space  coordinates 
of  electron   1.   In  particular,  the  one-electron  integrals 
programmed  and  calculated  in  the  course  of  the  present 
work  are  the: 
Overlap 

(1.4)  fx*   X'  dT  , 
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kinetic  energy 

1.5)  f  '^*^-   I  A)  X  '  dT  , 

nuclear  attraction 


(1.6)  J    %*   r;l  %'  dx  , 
dlpole  moment 

(1.7)  J  :i/  ""  r^  cos  9^^  X  '  dT  , 
quadruple  moment 

(1.8)  J   X*[r^(3  cos^  9^-   l)/2]  X'  dx  , 
electric  field, 

(1-9)  f  X""  v-^   cos  e^   X'  dx  , 

and  electric  field  gradient 

(1»10)  ft*   ^-^(3  cos^  Q^-   1)7l'  dx 


Integrals,  where  X  and  X  are  the  STO ' s  (II  4.2). 

The  only  two-electron  property  considered  here  is  the 
inter-electronic  potential   r~    in  the  electronic  energy 
equation.   For  this,  one  needs  the  integrals 
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(1.11) 


)*(1)  (t)j(2)  r'2  <t)i^(l)  <t^(2)  dT^  dTg 


5     C  .C  .C  ,  C  „  G 
^^^     pi  qj  rk  si   pqrs 


where 

''•'"'    Ws=ir^p'l'^q'2'  ^12  -^^rfD^s'^'  d-1^-2  ' 

The  Integrals  (1.12)  are  usually  divided  into  three 
types  according  to  the  centers  of  the  STO ' s  involved. 
This  division  is  relevant  only  when  considering  the  means 
of  calculating  or  tabulating  them.   The  three  types  are  the: 

Coulomb 

Hybrid 
(l.lt)  ff  Xlil)t^"'(2}    r-1  t;'(l)X;"(2)  dT,  dx,  , 

and  Exchange 


1.15)    /7X;(1)X;'(2)  r-l  X;'(1)X;"(2)  dx,  dxg 


integrals,  where  the  subscripts   a  and   b   denote  the 
centers   A  and  B  of  the  STO ' s .   The  integrals  having 
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all  AO ' s  on  the  same  center  can  be  obtained  from  any  of 
the  above  by  setting  R,  the  internuclear  distance,  equal 
to  zero. 

Most  of  the  literature  on  the  calculation  of  these 
integrals  and  the  earlier  calculations  for  some  closed- 
shell  systems  have  treated  the  real  STO ' s  which  are  defined 
as  the  real  and  imaginary  parts  of  the  STO ' s  (11  4.1). 
The  one-electron  integrals  are  the  same,  whether  one's 
notation  denotes  complex  or  real  STO's,  and  the  two- 
electron  integrals  of  complex  orbitals  are,  in  general, 
linear  combinations  of  the  integrals  of  the  real  STO's. 
Therefore,  the  STO's  discussed  in  this  section  and  in  the 
tables  referred  to  are  the  real  STO's. 

In  the  Coulomb  and  Hybrid  integrals,  the  charge 

distribution  for  electron   1,  X  ( 1 )  A-  ( 1 ) »  can  be 

a     a 

expanded  in  a  finite  series  of  real  spherical  harmonics 
on  center  A.   Then,  one  can  use  the  expansion 

I 

(1-16)   r-l  =  ri  ^1^  -^  HZ  S,^(0^^,  k{)    S^JG^^,    ^^) 
£=0  r^    m=-& 

where  the  S„   's  are  the  real  spherical  harmonics  and   r. 

and   r^   are  the  smaller  and  larger,  respectively,  of   r  ., 

and   r  ^.   The  integration  over  the  coordinates  of 
a2  ^ 

electron   1   can  be  carried  out,  yielding  a  finite 
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expression  for  the  potential  at  electron   2  due  to  the 
charge  distribution  of  electron   1.   The  integration  over 
(Pp   can  also  be  performed,  giving  a  finite  series  of 
integrals  over  the  two  remaining  coordinates  of  electron   2. 
In  the  next  two  sections,  two  alternative  means  of 
expressing  and  evaluating  the  integrals  appearing  in  this 
finite  series  will  be  given. 

In  the  exchange  integrals,  each  electron  has  a  charge 
distribution  which  is  a  product  of  an  AO  on  A   and  an  AO 
on  B   and,  in  general,  this  gives  rise  to  an  infinte 
series  of  two-dimensional  Integrals.   While  work  on  the 
Integrals  reported  here  was  in  progress,  other  workers 
produced  and  distributed  a  machine  program  for  calculating 
the  exchange  integrals  [25].   Hence,  plans  for  calculating 
and  Including  them  here  were  abandoned  and,  instead, 
programs  were  written  for  calculating  some  of  the 
properties  for  which  the  integrals  were  not  available. 

Three  sets  of  tables,  hereafter  referred  to  as  T  I, 
T  II,  and  T  III,  were  produced  as  part  of  the  present 
work  and  published  separately  [26].   Detailed  descriptions 
and  instructions  accompany  the  tables  and  will  not  be 
repeated  here. 
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2.   The   C^E^-Functions 

op 

The  spherical  polar  coordinates  used  for  defining 
the  AG's  {II  4.1)  are  as  described  in  Figure  1. 


Fig.  1  -Spherical  polar  coordinates  of  electron  at  point  P. 


All  of  the  integrals  (1.4)  -  (1.9),  (1.13),  (I.l4), 
and  most  of  those  in  (1.10)  can  be  expressed  in  terms  of 
the  integrals 


c^p  lP,t)  =  n 


oo 

K  r    2 


IT 


^a^^aj  ^^'^  ^a  ^^a  ^^P  ^  "^a^a  "  ^^b^b  ^ 


2.1 


a+7+£-l      B+5+e-l  7^  B^         •     e^         •    ^o 

r      ^  rr,  cos    0      cos   0,     sm   0      sm   0, 

a  b  a  b  a  b 


where 
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K  =  a+P+7+5  +  2£  +  l 


(2.2)  ^ 


t  =  (h^-h^)/(m.^  +  h^ 


The  integrals  (2.1)  exist  only  for  a  +  7+2£+1^0.   A 
change  to  elliptic  coordinates 


(2.3) 


4  =  (r  +  r,  )/R 
^     a   b  ^ 


n  =  (r^-^b^/R 


puts  the  integrals  (2.1)  in  the  form 

Cjp^(p,t)  =  (|)  J  d^^  dTi  exp  (-p4-pt  Ti) 
(2.4)  ^      -^ 

(?+Tl)''(?-Ti)^(l  +  |ri)^(l-  eTl)^(^^-  1)^(1-  Ti^)^ 


These  are   (l-t)~   times  the   C^g   -functions  defined  by 
Rildenberg  ,  Roothaan,  and  Jaunzemis  (hereafter  referred  to 
as  RRJ )  [27].   The  present  definition  is  formulated  so 
that  the  relation 

■y5£ 

can  be  used  to  simplify  tabulation  of  the   C^o   -functions, 
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RRJ  have  eiven   recursive  relations  for  the   C  □  - 
°  ap 

functions  which,  when  applied  to  numerical  values  of  the 
functions,  lose  accuracy  very  rapidly  due  to  "differencing 
effects".   In  a  later  paper,  [28],  Roothaan  applied  the 

recursive  relations  to  analytic  expressions,  of  a  simple 

■yBs 
general  form,  for  the   C''o   -functions  and  gave  formulas 

yB  £ 
for  all  of  the   Co   -functions  having  a  >  0  which  are 

ap  ^    — 

needed  for  the  one-  and  two-electron  Integrals  of  the 
Is,  2s,  2pcr,  and  2pTr  AG's.   A  single  subroutine,  requiring 
stored  arrays  of  polynomial  coefficients  for  each  C'g 
function,  was  written  and  used  successfully  to  compute  all 
necessary  C'^o   -functions  with  a  >  0.   The  procedure  is 
extremely  fast  compared  with  others  which  use  numerical 
integration  and  it  gives  high  accuracy  for  the  entire 
range  of  the  parameters.   By  other  methods,  the  sizes  and 
ranges  of  tables   T  I,  T  II,  and  T  III  would  have  made  the 
time  and  cost  of  their  calculation  prohibitive. 

yS  £ 

The  appearance  of  the  Cq      -functions  with  a  <  0 
^^  ap 

did  present  some  difficulties  since  they  are  of  a  more 
complicated  form,  involving  exponential  integrals  and  terms 
having  singularities.  In  order  to  avoid  these  complications 
in  the  calculation  of  the  Coulomb  Integrals,  a  method  was 

•y5  £ 

devised  for  expressing  the  linear  combinations  of   Co 

functions  appearing  in  them  with   a  <  0   in  terms  of 

y  6  £ 
C^Q   -functions  with  a  >  0. 
ap  — 
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By  using   (  ^+ri )    as  one  of  the  factors  in  two 
Integrations  by  parts,  with  respect  to   ^   and  r] ,    of 
C^g  ,  and  by  adding  the  two  results,  one  obtains 


(2.6) 


(a+1)  C^60  _^     750 
aP    [X    a+1, 1 


1  y     7-1,5,0     7,5-1,0  1  y,7+l,^ 
p  (^"  ^^a+2,P    +  ^^a+2,P   j"  ^aP 


The  functions  h  g   are  Integrals  having  the  property 


(2.7)  h^^  =  h^+2m,6+2n 

ap    a-2m,P-2n 


for  all   m,n   for  which  the  corresponding  integrals  in 
(2.6)  exist.   Then,  by  doing  the  same  thing  with   (^-ri) 

used  as  a  factor  in  the  partial  integration,  one  obtains 

■y  5 
a  formula  for  h  ^  , 

aP  ' 


(2.8) 


h^'^-^1  =  -(P+1)  c^^°  -^c^'o 
aP        "^      aP    |j.   a,P+l 


1  r   7-1,5,0     7,5-1,0  1 
+  p  [^^a,P+2    -  ^^a,3+2   j 


Equation  (2.6)  gives  the  linear  combination  of   C'^g 
functions  on  the  left  side  of  (2.6)  in  terms  of  those  with 
higher  a.   This,  with  formulas  for  raising   e,  can  be  used 

75  E 

to  convert  the  linear  combinations  of   C'q   -functions  for 

ap 
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the  Coulomb  Integrals  Into  expressions  involving  c'^a^    's 

with  a  ^  0   and   ho  's  with  a  <  0.   Then,  (2.7)  can  be 

used  with  convenient  choices  of   m  and   n   to  obtain  the 

ho  ' s  In  terms  of  h^o  's  with  a  >  0,   The  latter  are  then 
ap  ap  — 

given  by  (2.8)  In  terms  of  C^ a        with  a  >  0.   Relations 

ap  — 

between  the   C  g   -functions   can  then  be  used  to  simplify 
the  final  results. 

•y5  £ 

The   C  Q   -functions  with  a  <  0  which  are  needed  for 
ap 

some  of  the  one-electron  and  Hybrid  Integrals  were  calcula- 
ted rapidly  by  using  explicit  expressions  given  by  Roothaan 
[29],  but  these  expressions  are  subject  to  serious 
differencing  errors  In  the  range  of  parameters  corresponding 
to  small  pt   (see  Equations  (2.2)).   Therefore,  for  these 
ranges,  and  for  some  of  the  Integrals  (1.10)  which  cannot 

y  5  £ 

be  expressed  In  terms  of  Cq      -functions,  the  slower  but 
^  ap  ' 

more  accurate  methods  of  the  next  section  had  to  be  used. 


3-   Barnett-Coulson  Zeta-f unctions 

The  basis  of  the  Barnett-Coulson  Zeta-functlon  method 
[30]  is  the  use  of  the  expansion 

-R  _i 

r^-1  e^  ^^  =  p-+l^  (2n+l)(tT)"2 

(3.1)  ""^^ 

P  (cos  0  )  C   (l,t;T) 
n      a   ^mn   '  ' 


See  Page  213,  Ref .  [27] . 
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where 

t  =  Pr^ 

T  =  PR 


and  where   P   is  the   n-th  Legendre  polynomial.   The 

functions   C    are  defined  by  recursive  relations  which 
^mn  '^ 


start  with  the  functions 


3.2)  ^On  =  ^    1  (tj  K    ^  (tj 


n  +  2       n  +^ 

where   t^   and   t^   are  the  larger  and  smaller,  respectively, 

of   t  and   T  and  where   I    -,   and  K    -,   are  the  Bessel 

n  +  "2        n  +  "2 

functions  of  purely  Imaginary  argument-   The   C„j^   functions 
with  m  =  0,  1,  and   2   are  denoted  by  7^,  p^,  and   q^. 

Using  (3.1)  and  trigonometric  Identities,  the  co- 
ordinates of  an  AO  on  center  B   can  be  expressed  in  terms 
of  spherical  polar  coordinates  on  A  which  may  then  be 
used  as  the  variables  of  integration.   This  enables  one  to 
express  all  one-electron  Integrals  and  the  two-electron 
Coulomb  and  Hybrid  integrals  in  terms  of  the  one-dimensional 
integrals  defined  by  Barnett  and  Coulson, 

00  ^  _^1 

(3.3)     Z       iC^t)  =  f    e-'^^  C^n(l,t;T)  t    ^  dt  , 
m,n,^  +  2        "q 
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where  K  and  x  depend  upon  the  orbital  exponents  and  R. 

The  Zeta-functions  (3-3)  with  m  =  0,  1,  and   2   are 

denoted  by  G     ,  ,   P     -,  ,  and   Q     -,  respectively. 
n,i+-2     n,^  +  -p         n,i+-2 

In  what  follows,  functions  are  defined  and  formulas 
are  written  for  the  express  purpose  of  using  them  for 
calculation  with  an  automatic  computer.   Their  formulation 
is,  therefore,  determined  so  as  to  minimize  the  number  of 
machine  operations,  avoid  overflow,  and  yield  maximum 
accuracy.   Consequently,  some  functions  and  formulas  do 
not  appear  in  the  most  familiar  or  convenient  way  of 
expressing  them  in  ordinary  mathematical  text. 

The  functions 

1    _   1 

(3.^)      B^(z)  =  7r^(2z)     ^(2n+l),'(nJ  )"^  I    ^(z)e"^ 

n+  "2 

--     n  +  - 

(3.5)  C^(z)  -  TT  ^(2z)    2[(2n+l).']"^  n,'  K    -^(z)  e^ 

n  +  -2 

i  ^>-^< 

(3.6)  ^mn^^^''^  =  ^^^^'      ^rm^^'^'''^    ^ 


will  be  introduced  and  7    ,    P    >    and   q   will  be  used  to 

^n   ^n       ^n 

denote   C    for  m  =  0,  1,  and   2,  respectively.   By 
^mn  '   '       J     ^ 

these  definitions, 

(3.7)  7^{t,T)   =    {t/T)^{tJtS     ^  B^^tJ    C^(tJ 
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and  Barnett  and  Coulson's  Zeta-f unctions  (3.3)  are  just 

-1/2 
T      times  the  functions 

(3.8)      Z^^^C^t)  =     J      e  C^n^t'-^)  ^   <^t 

0 

with  the  same  values  of  m,  n,  and  £.      The  z"^ „      functions 
with  m  =  0,  1,  and   2  will  be  denoted  by  G   ^ ,    ?    . ,    and 
Qy^n,    respectively.   The  recurrence  relations  for  the  barred 
functions  defined  here  are  the  same  as  those  given  by 
Barnett  and  Coulson  for  the  corresponding  unbarred  functions, 
The  first  of  the  recurrence  relations  for  raising  m  is 


(^•9)      ^ni   =  -(2n+l)-l  (Gn-l,i+l  "  ^^+1,^+1) 


For  convenience  in  deriving  and  expressing  integrals 
of  the  AO ' s  the  functions 


,  aVn.+£+l       p       o  r  -ctr    -  Pr, 

z"^(K,t)    .  ^-3_      /  r^dr  /sine      d0      e        ^  ^ 

nz      '  2  ^/         a     a  J  a        a 

0  0 


(3.10) 


r  r,         cos      9      P      cos   6        , 

a  b  an  a      ' 


where   K  =   ct/P,  will  be  introduced  here.   For  k  =   0, 
these  are  equal  to  the  functions  (3.8).   A  familiar  relation 
among  Legendre  polynomials  gives  the  formula  for  raising 
the  index  k, 
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— mk     I  r,    ..N-li"  77m, k-1  ,  /  , -,  sTT-m.k-l  1 
(3.11)     Z^^   =   (2n+l)   (nZ^^^^^  +  (n+l)Z^;^^^  j-  . 

By  usine;  one  of  the  recurrence  relations   for  the   C 
functions,  one  can  obtain  the  formula  for  raising  m. 


-^  nH  ni       n,^+2       n,i+l 


To  calculate  the  Zeta-f unctions  {ji.Q) ,    the  recurrence 
relations  (3' 9)  and  (3-12)  in  m  are  used,  starting  with 
an  array  of  the  G  .   functions.   This  produces  no 

X\xj 

appreciable  loss  in  accuracy  for  the  integrals  considered 
here.   The  computation  of 

r"  -'<t-{t,  - 1,)  _  0 

(3.13)     ^ni^'^'^^   "  J   ^  7n(t,T)  t*  dt 

0 

is  performed  by  dividing  the  range  of  integration  at   t  =  t 
and  expressing  the  result  in  the  form 

(3.1^)     G^^C^'-^)   =   ^n^""^  ^ni^'^'""^  +  ^n^''^  ^ni^'^'''^  • 
By  a  change  in  the  variable  of  integration,  one  gets 


(3.15)     G%   =   T^+1   /   e-^t^'^-l^^+l]  B  (Tt)  t^^^-^1  dt 

0 


See  Page  224,  Ref.  [3O] 
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oo 


(3.16)    G^^  =   x^+1  f    ^-T[i^c+l)t-l]    ,^(^^)  ^i-n^^  ^ 


1 

."   ill   UllC     V 

functlons,  are  calculated  by  the  formulas 


The  B   functions,  which  also  appear  in  the   C' □ 


B_,  =  e~   cosh  z 
(3»17)         B„   =  z~   e~   sinh  z 


B   =  (4n^-l)  z"^  (B   „  -B   J  . 
n  n-2   n-1 


For  low  values  of   z,  the  forward  recurrence  procedure 
loses  accuracy.   Therefore,  for  z   below  a  certain 

criterion,  the  backward  recurrence  procedure,  starting  at 

t     I 
high  n  with  B  =  B   ,   equal  to  an  arbitrary  number, 

is  applied,  giving  a  series  B  , -,  ^  •  •  •  5  B  ,  .   Then  the 


correct  values  of  the  B   functions  are  obtained  by 

n 

multiplying  the  whole  series  by  B  -|/B_-,   where   B  .,   is 
calculated,  without  differencing  errors,  from  the  first 
equation  in  {J).17). 

The  recurrence  scheme  for  the   C    functions  is 

n 


C  ^  =  z  ^ 


{3.18)  ^0  =  ^ 

_^  =  (2n+3)~"^  {(2n+l)C^  +  z^(  2n+l )  ~^C^_^} 


n+J 
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This  procedure  is  accurate  for  Increasing  n  and  all  z  _^  0, 

The  G"""  „   functions  lose  accuracy  rapidly  when 
nJi 

recurrence  relations  in  n  and   i   are  applied,  so  these 
are  calculated  by  Legendre-Gauss  numerical  integration. 

For  the  G^.   functions,  recurrence  methods  work  quite 

n  i/ 

well,  so  the  following  procedure  was  devised  for  their 
calculations : 

First,  the  functions  with  n  =  0   are 

00 

(3-19)    Gq^  =   e   t     /   e         t  dt  . 


By  partial  integration,  the  recurrence  relation 


-s      ,,,^,  ^-l  5  £      -Kt 


(3.20)     Gq,   =   (K+1)"^  It      e""''  +  ^G 


is  obtained.   This  and  the  starting  function 

(3.21)  G^Q   =   (K+1)"^  e"'^'' 

give  Gq^      for  all  £  ^  0.      For  £    <   0,  (3-20)  is  used  to 
decrease  £,    starting  with 

(3.22)  G^^_^  =  e-"^  f-e-^^'^+l)  Ei(-T(K+1))} 
where 

00 

(3-23)     Ei(z)   =  -  /   t~^  e^'^  dt 
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— g 

is  the  familiar  exponential  Integral.   By  the  form  of   G  ^ 
in  (3.16)  and  the  definition  of   C_-,   and   C„   in  (5.18), 
it  is  seen  that 

(3.24)       a\^^    ,  G=^^  T-i  . 

The  above  procedure  provides  all  the  necessary  G^.  's 
with  n  =  -1   and   0,   For  higher  values  of   n,  the 
relation 

(3.25)   gJ^  =  (2n+l)-l  {T(2n-1)G^_^^^_^  +  x^l 2n-l ) -1g^_2^ ^ 

is  used.   Equation  (3.25),  which  is  derived  by  substituting 
(3.18)  in  (3.16),  does  not  lose  accuracy  as   n   is  increased. 


4.   Tables  of  Integrals 

Tables  T  1,  T  II,  T  III  can  be  used  to  obtain  numerical 
values  of  the  one-electron  integrals  (1.4)  -  (I.9),  most, 
but  not  all,  of  the  electric  field  gradient  integrals 
(1.10),  and  the  two-electron  Coulomb  (I.13)  and  Hybrid 
(l.l4)  integrals.   The  range  of  parameters  is  such  that  one 
can,  by  interpolation,  obtain  approximately  six  significant 
figures  in  the  integrals  required  for  the  calculation  of 
the  first  row  diatomic  molecules  for  a  range  of  R  at 
least  up  to  their  dissociation  distances. 
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The  Integrals  in  T  I   can  be  used  with  the  simple 
formulas  given  there  to  obtain  the  overlap,  kinetic  energy, 
and  nuclear  attraction  integrals.   With  the  exception 

already  mentioned  above,  the  integrals  (1.4)  -  (1.10)  can 

■y  5  £ 
also  be  obtained  in  terms  of  the  C   a        functions  of  T  III. 

ap 

Explicit  expressions  for  these  and  a  number  of  other  one- 

■yb  £ 
electron  integrals  in  terms  of  the  C  q   functions  can 

°  ap 

easily  be  written  by  referring  to  the  definition  (2.1)  of 

the   C^Q    functions, 
ap 

Aside  from  a  simple  multiplicative  factor,  two-electron 
Coulomb  integrals  of  the  real  STO ' s  are  given  directly  in 
T  II   as  functions  of  p   and   t   (see  Equation  (2.2)).   All 
Coulomb  integrals  for  the  real  or  complex  STO ' s  defined 
above  can  be  expressed  as  simple  linear  combinations  of  the 

integrals  in  the  tables.   To  calculate  these  integrals, 

y  b  £ 
formulas  for  them  in  terms  of   C^q    functions  with  a  =*  0 

aP  — 

were  derived  by  the  method  described  in  Sections  1  and  2. 
In  deriving  formulas  for  the  Coulomb  integrals,  one  can 
perform  the  first  integration  over  the  coordinates  of  either 
electron  1  or  2.   Thus,  one  can  form  two  different  expressions 
for  each  non-symmetric  Coulomb  integral.   The  values  of  the 
argument   t   in  the  two  formulas  are  opposite  in  sign  and, 
in  general,  the  expression  having   t  ^  0  is  subject  to 
differencing  effects.   Therefore,  two  expressions  for  each 
Coulomb  integral  were  programmed  and,  in  all  cases,  the  one 
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with   t  ^  0  was  used.   The  tabulation  of  the  Hybrid 

Integrals  Is  more  complicated  since  they  are  functions  of 

three  parameters  even  after  a  multiplicative  factor  is  taken 

out.   Therefore,  tables  of  the  two-parameter  C  g   functions 

in  T  III   were  constructed  for  their  calculation.   The 

user  is  then  required  to  obtain  appropriate   C  g   functions 

from  the  table  and  form  linear  combinations  of  them  for  each 

Hybrid  integral. 

As  mentioned  previously,  the  analytic  expressions  for 

C'^g^,  ct  <  0,  provided  by  Roothaan,  do  not  give  accurate 

results  in  the  range  of  parameters  where  pt   is  small. 

Therefore,  this  procedure  was  used  for  pt   large,  where 

accurate  results  could  be  obtained,  and  the  slower  Zeta- 

function  method  (Section  3  above)  was  used  to  calculate  the 

C'^'q^  's  for  the  rest  of  the  table.   Fortunately,  the  latter 
ap 

method  is  faster  and  more  accurate  in  precisely  the  regions 
of  the  table  where  the  former  gives  errors.   There  were 

•y  5  £ 

large  regions  of  each  table  of   C''g   functions  where  both 
methods  gave  the  same  results. 

After  the  tables  were  finished,  the  programs  used  in 
their  calculation  were  converted  into  subroutines  which  are 
capable  of  producing  all  of  the  integrals  (1.4)  -  (1.10). 
For  a  given  set  of  parameters,  these  integrals  can  be  cal- 
culated in  a  few  seconds.   The  programs  for  the  calculations 
reported  below  used  these  subroutines  for  all  the  integrals 
involved. 

-87  - 


IV  A  NWIERICAL  METHOD  FOR  SOLVING  THE  RADIAL  3CHR0EDINGER 
EQUATION 

1,   Introduction 

The  wave  equation  for  the  nuclear  motion  of  a  diatomic 
molecule,  In  the  Born-Oppenheimer  approximation.  Is  one 
which  Is  encountered  frequently  In  quantum  theoretical 
calculations.   Numerical  methods  for  Its  solution  have  been 
developed  and  used  [L|-, 31, 32,33  J  over  many  years  for  atomic 
problems  where  the  potential  is  one  obtained  by  Hartree-Fock 
self-consistent  fields  or  the  Thomas-Ferml-Dirac  statistical 
field  methodso   Only  relatively  recently  have  computational 
techniques  and  the  application  of  electronic  computers  enabled 
one  to  obtain  accurate  theoretical  internuclear  potentials  at 
enough  Internuclear  distances  to  calculate  the  wave  functions 
for  the  motion  of  the  nuclei  and  use  them  to  obtain  averages, 
over  nuclear  motion,  of  molecular  propertieso 

The  present  investigation  is  concerned  with  obtaining 
an  accurate  method  for  calculating  the  nuclear  wave  functions 
and  vlbrational-rotatlonal  energies  of  diatomic  molecules 
with  some  economy  in  the  number  of  values  of  the  inter- 
nuclear potential  required.   In  this  chapter,  an  improved 
formula  for  the  correction  of  trial  eigenvalues,  which  does 
not  depend  so  much  for  its  accuracy  upon  the  smallness  of  the 
step-size  in  the  radial  coordinate,  and  an  analysis  of  the 
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converpence  of  tb8  nrccedure  are  given.   A  computer  sub- 
routine was  written  and  numerical  results  obtained  froin 
it  are  fiven  for  a  cpse  where  exact  analytic    solutions 
are  known. 

In  what  follows,  the  vibrational  quantum  number  v, 
V  =  0,1,2,...^  will  be  used  as  a  subscript  to  index  the 
eip-envalues   S   v;ith  the  usual  convention  that 

V 

3   <  E,  <  Ji'  <  ,  ,  .  . 
o  —  1  —  ^  — 

?. .   Method  of  InteFrstion . 

The  ochroe dinger  wave  equation  for  the  motion  of  the 
nuclei  of  a  diatomic  nclecule.  rep-srdea  as  a  symmetric  top, 
can  be  expressed  in  polar  spherical  coorainatos   i^j®,!: 
of  one  nucleus  relative  to  the  other  and  one  additional 
coordinate,   (j),   pfivinp  the  angular  orientation  of  the 
electronic  charge  cloud  about  the  internuclear  axis.   It 
has  been  shown  L3il-J  that  the  wave  equation  can  be  separated 
into  its  ane-ulai'  and  radial  narts  and  its  solution  may  be 
expressed  in  the  form 

where  ?{Y{]      is  tne  solution  of  the  one-dimensional  radial 
Schroedinger  equation 

(2.1)  F^^^R)   =   (II(R)  -5)  P(R)  , 

P^^^R)   =   d''p(R)/dR'', 
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the  Y^. „'s   are  the  hyper geometric  functions,  A  is  the 
z-component  of  electronic  angular  moinentum,  and  the  radial 
internuclear  potential,   U(R),   is  of  the  form 

u(R)  =  [j(j+i)  -^A?]  r"^  +  z  z^r"-"-  +  E  ,(R)  . 

a  D       ei 

The  second  term  is  the  electrostatic  Coulomb  repulsion 
energy  of  the  nuclei  and  E  -,  (R)   is  the  electronic  energy 
obtained  by  solving  the  electronic  wave  equation  discussed 
in  Chapter  II.   The  boundary  conditions  for  (2.1)  are: 
P(0)  =  0,   P(R)   bounded. 

To  get  a  difference  equation  whose  solutions  approxi- 
mate the  solutions  of  (2.1)  we  let 

R^   =   Ih     ,         1  =  0,1,2,  .-,n+l, 
(2.2)  P^   =  P(R^)   , 

U^   =   U(R^)   . 

By  dropping  terms  of  the  fourth  order  and  higher  in  the 
series 

and  using  the  differential  equation  to  replace   P.    ,   one 
obtains  the  simple  integration  formula 

(2.1,)      P^^^  +  P._^  -  2P.   =   h2(U.-E)P. 

which  has  an  error  of  approximately  yp  -^-i   "  ^   higher 
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order  Integration  formula,  which  does  not  Involve  any  more 
values  of   P.,   can  be  obtained  by  subtracting  h  /l2   time 


the  series 


00  o.2k    (2i,+2) 


from  {2.3) o      Then,  dropping  sixth  and  higher  order  terms 
in  h  gives 

(2.6)      Y,.T  +  Y,  ,  -  2Y.   =  h^(U.  -E)P. 
1+1     1=1      1         1      1 


where 


2 


(2o7)   Y^   =  Pi  -X2  ^1^^   "   [1-^  (U^-  E)]P^   c 

The  error  in  (2 ,,6)  is  approximately   -  'oV^   P-s   ■>   "^^^ 
integration  formula  (2,6),  usually  attributed  to 
Numerov  [35],  involves  the  extra  calculation  of  the  Y.'s 
but  reduces  the  number  of  points  where  values  of  U(R) 
are  required. 

In  the  range   E  >  U(oo),   corresponding  to  unbound 
states  of  the  two  nuclei,  solutions  of  (2ol)  exist  for 
all   E   and  can  be  approximated  by  simply  using  (2o[|.)  or 
(2o6)  to  integrate  outward,  starting  with  the  boundary 
values. 


(2.8)     P   =0,     ^1~^  small  arbitrary 


number. 


For  E  <  U(oo),   the  two  nuclei  are  bound  together 
and  solutions  of  (2.1)  exist  only  for  a  set  of  discrete 
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values  of  E,   the  eigenvalues  of  the  problem.   The  method 
for  calculating  these  eigenvalues  and  corresponding  solutions, 
or  eigenfunctions,  is  the  subject  of  the  present  chapter. 

The  boundary  condition,   P(R)   bounded,  is  approximated, 
as  usual,  by  the  conditions 

P  ^^   =  a  small  arbitrary  number 

(2.9) 

P   =  P  ,  T  exp(R^,T  yu  ,,  "E  -  R  ■/U~^)    , 
n      n+1   ^  n+1  "    n+1      n  ^    n 

The  second  of  these  conditions  results  from  the  assumption 

that,  at   R  ,   U(R)   is  slowly  approaching  a  constant. 

The  usual  numerical  procedure  is  to  provide  a  first 

estimate  of   E   and  integrate  outward  from  R  =  0   to  some 

point   R   using  starting  values  (2.8)  and  integration 

formula  (2,14.)  or  (2,6).   Since  (2.1)  is  homogeneous  in 

P(R),   the  resulting  P.   values  may  be  replaced  by 

P^^*  =  Pj^/Pj^,   1  =  l,2,,..,m.   With  starting  values  (2.9), 

the  same  procedure  is  used  to  integrate  inward  from  R  ^ 

°  n-1 


m 


ch 


to  R        and  yields  values   P.  ,   i  =  n-l,n, ...,m   su 

that   P-^^  =  P°^*  =  lo   Then,  a  correction  to   E  is  deter- 
m     m  * 

mined  by  the  difference  between  the  slopes  of  the  two 
curves  at  the  crossing-point   R   and  the  process  is 
repeated  until  the  two  curves  meet  with  the  same  deriva- 
tive.  The  correction  D(E)   is  usually   calculated  by 

See  Equation  6,  page  86,  of  Ref.  [k.] , 
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using  the  formula 

00 

(2.10)      D(E)   =   (Pout'^in^  //p(R)  dR 

where  the  terras  in  the  numerator  are  the  derivatives  at 

R    of  the  curves  resulting  froin  the  outward  and  inward 
m  ° 

integration,  respectively. 

Equation  (2.10)  has  been  derived  from  the  differential 
equation  and  is  claimed  to  give  first  order  convergence  in 
the  error.  The  means  that  if  a  trial  solution  satisfies  the 
differential  equation  (2.1)  at  all  points  except   R  ,  then 

(2.10)  deviates  from  the  true  correction  by  an  amount  which 
is  proportional  to  the  correction.   In  the  next  section,  it 
will  be  shown  that  the  eigenvalues  of  the  difference  equa- 
tions, (2.I4.)  or  (2.6)  are  the  zeroes  of  appropriate  functions 
of  the  trial  eigenvalue   E  and  that  these  zeroes  can 
conveniently  be  calculated  by  the  Newton-Raphson  method. 

The  correction  is 

(2.11)  D(E)   =   -  P(E)/P' (E) 

where   P(E)   is  a  function  whose  form  depends  upon  the 
integration  formula  used.   For   E  near  an  eigenvalue  E  , 

(2.12)  D(E)   =  ^E  +  (/\E)^Q  +  higher  order  terms  in  /_^E 

where   Ae  =  E  -E   and 
^— ^     o 

(2oi3)         Q  =  f"(e;^/2p'(e;;  . 

-  93  - 


Thus,  the  convergence  of  the  present  method  is  of  second 
order  in  /\E.   In  Section  L|.,  an  analysis  of  the  behavior 
of   D(E)   over  the  whole  range  of   E   is  given. 

3.   Correction  Formula. 

The  integration  formulas  (2.1j.)  and  (2,6)  can  be 
written  as  systems  of  equations  in  the   ^^'s   ^-^^  Y  's 
respectively,  as  follows: 

(2h"^  +  U^-  E)P-L  -  h'^Pg  =  0 

(3.1)   -  h"^P^_-|_  +  (2h'^ +U^ -E)P^  -  h"^P^^-,_   =   0 


and 


=  0 


{2h"2+(U^-E)[l-  I2  (U3_-E)]'^?Y^-h2Y2 

2 
(3.2)  -h"^Y^_^+f2h"^+(U^-E)[l-^  (U^-E)]"^]Yj^-h"^Y^^3_  =  0 

2 

-h'^„  ,+  fh"^+(U  -E)[l-^  (U  -E)]"^}y„         =  0 
n-1   ^      n       Id       n      ^    n 

where   i  =  2,3,».«,n-l   in  both  cases.   A  term  involving 
E  has  been  omitted  from  the  n-th  equation  in  both  (3»1) 
and  (3«2).   It  can  easily  be  shown  that  if  one  has,  as 
one  should,  started  the  inward  integration  at  a  large 
enough  R   to  justify  the  assumption  used  in  forming  the 
boundary  condition  (2.9),  this  omission  is  justified. 
The  procedure  used  here  to  derive  and  analyze  the 
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correction  formula  is  an  adaptation  of  a  very  general  and 
effective  technique,  due  to  LBwdin  [36J,  for  calculating 
solutions  of  the  Schroedinger  equation.   To  apply  L8wdin's 
method,  consider  the  vector-matrix  fonnulation  of 
equations  (3.1)  or  (3«2), 

i3>3)  M  C  =  £   , 

where   0   is  a  null-vector,   _C   is  a  vector  containing  the 
P  's   or  the  Y  's,   and  M   is  the  symmetric  matrix  of 
coefficients  which,  for  the  present,  may  simply  be  regarded 
as  functions  of   £.   After  an  outward  and  an  inward  integra- 
tion with  a  trial  value  of   E,   all  equations  except  the 
m-th  are  satisfied.   Now,  assume  that  the  first  row  of  M 
and  the  first  element  of   C_  correspond  to  the  m-th 
equation  and  the  m-th  variable,  respectively,  in  equations 
(3.1)  or  (3»2).   Then,  by  partitioning  the  first  index  of 
M  and  C,      the  result  of  the  integration  can  be  expressed 

^11   ^Il^/l  \     /f(e)\ 
(3.1+) 


^al   ^aa 


V^  / 


where      P(E)      is    the   amount   by  which  the  m-th   equation  of 

(3.1)    or    {3'2)    is   not   satisfied  when   integrating  with  the 

trial    value      E.       It    is    assumed   that    the   first    element    of 

C,      which  is      P        or     Y    ,       is   non-zero   and,    for    convenience, 
—'mm'  '  ' 
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it   is    set   equal   to   1.      Equation    (3«i;)  niay  be  written  in 
the   form 

(3.5)  P(E)      =      M,,    +  Mtn    C^ 

11        -al    -a 

(3.6)  0     =     M„,    +  M„„   C„    . 

—  -al        -aa   -a 

The  function  P(E),   whose  zeros  are  the  eigenvalues  of 
the  difference  equations,  is  defined  by  (3»5)  in  terms 
of   C   which,  in  turn,  is  defined  as  a  solution  of  (3.6). 
Furthermore,  if  E   is  such  that   |M   |  7^  0,   then  the 
solution  of  (3 '6)  is  unique  and  P(E)   is  uniquely  defined, 
An  expression  for   F  (E)   is  obtained  by  differentiating 
(3.5)  and  (3.6)  with  respect  to   E   and  obtaining 

(3.7)  f'(E)   =  m;.^  +m;^  5a  ■'^l  ^-l 

Multiplying  (3.8)  on  the  left  by   c"*"   gives 

(3.9)  0   =   C"^  m't  +  C"^  m'   C„  +  C"*"  M„„  c'  . 

—     -a  -al   -a  -aa  -a   -a  -aa  -a 

Prom  (3.6)  it  is  seen  that  if  (3.?)  and  (3.9)  are  added, 
the  terms  containing   C    cancel  and  the  result  is 

(3.10)  P'(E)   =  m'   +  M't  C   +  C"*"  m't  +  C"^  m'   C^ 

11   -al  -a   -a  -al    -a  -aa  -a 

which  can  be  written 


-  96  - 


(3.11)  f'(e)    =    c"^  m'  c 


In   the    two    cases    considered   here,    only   diagonal 
elements    of     M      depend   upon     £      so    that   the    correction 
formula    (2.11)    can  be   written 

(3.12)  D(E)      =      -   FiE)/YZ  C?   m!.       . 

1=1      ^      ^^ 

In  the  case  of  (3.1),  the  derivatives  of  the  diagonal 
elements  of  M   are  all  equal  to   -1   so  (3«12)  becomes 

(3.13)  D(E)   =   [(-P^  ,  +  2P^-P^  n)h"^+(U  -E)P^J/2Z  P?  • 

m-i    m   m-±       m         m  Y=t"  i 

It  is  of  some  interest  to  compare  this  with  the  correction 

formula  (2.10).   If,  in  (2.10),  central  differences  at   R 

'  '  m 

are  used  to  estimate  the  derivatives  and  the  trapezoidal 
rule  is  used  to  estimate  the  integral  in  the  denominator, 
then  (2.10)  is  the  same  as  (3*13) •   Therefore,  use  of  the 
integration  formula  (3»l)  and  correction  formula  (2.10) 
as  described  above  yields  a  second  order  process  in  /\E» 

In  the  case  of  the  Numerov  integration  formula,  where 
equations  {3'2)    are  solved,  the  derivatives  of  the  diagonal 
elements  of  M   are 

2 

(3.11^)   m]_^   =   -[1-^2  (U^-E)]"^  ,      i  =  1,2, ..c,n. 

By  substituting   P . ' s   for   Y . ' s   the  correction  formula 
(3.12)  becomes 
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(3.15)  D(E)  = 


?  n   p 

[(-Y^  ,+2Y^-Y^  ,)h  +(U  -E)P^]/ZI  Pf 


h'      Convergence . 

A  rigorous  treatment  of  the  convergence  of  the  method, 
when  the  Numerov  integration  formula  {3'^)    Is  used,  would 
be  quite  complicated.   Furthermore,  since  the  difference 
in  the  two  methods  considered  is  a  term  of  order  h^  in 
the  solution  P(R)   and  since  both  are  second  order 
processes  in  /\E,   one  is  justified  in  assuming  that  the 
essential  features  of  the  convergence  of  the  two  procedures 
will  be  about  the  same.   Therefore,   the  convergence  of  the 
method  using  equations  (3.1)  will  be  considered  here. 

To  get  expressions  for   P(E),   f'(E),  and   D(E),   let 


(Il-.l)  H  = 


/_   2?i 

K. 

K 

•               . 

K         1 
K 

K 

"       ^ 

K 

K 

m-1 

K 

\ 

m+1    ■_ 

1 

K 

K 

G 
n 

where 
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K  =   -  h"^   , 

G.      =      2h"^  +  U.  ,       i  =  l,2,..,,n-l, 

G   =  h"^  +  U    , 
n  n   * 

and  where  the  blank  spaces  denote  zero  elements. 
Then,  for  equations  (3.I), 

(I4..2)  M  =   H  -  1  E 

where   1   is  the  identity  matrix.   The  symmetric  matrix 

M    can  be  diagonalized  by  an  orthonormal  matrix  U, 

whose  columns  are  the  eigenvectors  of  both  M„„   and  H„  . 

^  -aa       -aa 

Thus, 

(i|o3)              M    =  U  A  U"^ 
^  -^  -aa     

where  ^  is  a  diagonal  matrix  with  the  eigenvalues  of 

M    on  its  diagonal.   Then,  if  no  eigenvalue  of  M    is 
-aa  °  '         °  -aa 

zero,  a  solution  of  (3.6)  exists  and  can  be  written 

Letting 

(I4..5)  I  =   U"*"  %i 

and  substituting  (Li-.ij-)  for   C    in  (3*5),  one  obtains 
(1+.6)  P(E)   =   M^^  -  V"^  a"-'-  V  , 

11    '^ —   V   V 

V 
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As  one  can  see  from  il\..2) ,    the  diagonal  elements  of  A  are 


(1+.7) 


=   e   -  E 

V 


where   e  ,   v 
v' 


a-A^A9»*«  ij."*^  « 


e   <  e  J.T  .   are  the  eieen- 
V  —  v+1*  ^ 


values  of  H 


-aa' 


the  matrix  obtained  by  deleting  the  first 


row  and  column  of   H.   Hence, 


(1+.8)   P(E)   =   2h"^  +  U^  -  E  -  XI  vj  (e^-E)'^ 


and 


(1+.9)  f'(E)   =   -  1  -  Z  %  (^v"^^ 


-2 


The  set  of  points   R,  ,...,R   -,  ,   R  ,t,...,R  , 
^  1*         '    m-1'    m+1*    '  n' 

referred  to  by  the  subscript   "a"   in  the  vectors  and 

matrices  above,  may  be  partitioned  into  the  set  of  points 

Rt,...,R   -,  ,   used  in  the  outward  integration,  and  the  set 
1*    *  m-1'  ^ 

of  points  R  ^, ,o..,R  used  in  the  inward  integration. 
Letting  "out"  and  "in"  subscripts  denote  the  subvectors 
and  submatrices  resulting  from  this  partitioning,   H 

*"  9.0. 

can  be  written 


(I^.IO) 


H 


aa 


-out 


v2 


0 


-in 


where   H   .   and   H,    are  the  block  matrices  lying  on  the 
-out       -in  "^   ^ 

diagonal  in  (l+.l)  and  enclosed  by  the  dashed  lines.   By 
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partitioning  in  this  manner,  {L|..3)  can  be  written 


(l^..ll) 


M.=U.A,U. 
-out     -out  -out  -out 


M.    =  U.  A,       uT 
-in      -xn  -in  -m 


where  the  columns  of   U   .   and  U.    are  the  eigen- 

-out       -m  ^ 

vectors  of   H   .   and  H.    respectively.   Therefore, 
-out       -xn     J-       J  > 

e  ,   V  =  l,2,...,n-l,   are  the  eigenvalues  of  §   . 

and   H.  ,   Prom  the  form  of   H,   it  is  evident  that 
-xn  — ' 

these  are  the  solutions  which,  on  the  outward  and  inward 
integrations,  respectively,  lead  to  P   =  0  which  is 
contrary  to  the  assumption  that  the  first  element  of   C 
is  nonzero. 

To  show  that   P(E)   is  undefined  at  and  only  at  the 
eigenvalues   e    it  is  necessary  to  show  that  no   V   is 

zero  in  (ii.  8).   This  is  evident  since  ([|..5)  and  the  form 

-2 
of  M  ,   imply  that  the  elements  of   V  are   -h    times 

the  elements  in  the  last  row  of  U   .   and  the  first  row 

-out 

U.  ,   If  any  such  element  of   U   .   or  U.    were  zero, 
-xn        ''  -out      -xn  ' 

then,  since  M   .   and  M.    are  tridiagonal  matrices, 
'        -out       -xn  ^  ' 

the  column  containing  that  element  would  be  zero.   This, 
of  course,  is  not  so,  since  these  columns  are  the  eigen- 
vectors of  M   .   and  M.  . 
-out  -xn 

It  may  be  of  some  interest  to  show  how  the  character- 
istic polynomial, 

-  101  - 


(i;.12) 


P(E)   =   |M|   , 


of  (3.1)  is  related  to   P(E).   From  (L1..3)  and  (1;.5),  one 
gets 


(i|.13) 


P(E)   = 


Mil    ^ 


V 


A 


which  can  be  expanded  in  the  elements  of  the  first  row 
and  column  of  the  determinant  to  give 


{k.,lk) 


p(E)    =    M,^  T~r  ^  ,  -  T~  v^  TT    ^  ^ 

11  '  J   v'    -^ V  '  ,  V   V 

v'  V  v^fv 

From  (I4..6)  and  (I+.7), 


11   •' V   V  ■   '   '   V 

V  V 


(1^.15) 


P(E)   =  F(E)  TT  (e  -E)   . 

V 


Equations  (l+.S)  and  (I|..9)  enable  one  to  determine  the 
general  behavior  of   F(E)   and   D(E).   They  show,  first  of 
all,  that   F(E),   F  (E),   and   D(E)   are  defined  and 
continuous  for  all   E   except   E  =  e  ,   v  =  l,2,...,n-l 
and  that 
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(ij..l6) 


F(E) 

<iJ 

211"^   +   U      -   E 

m 

for 

E| 

large. 

P(E) 

< 

2h"^  +  U     -   E 

for 

E  < 

ei 

F(E) 

> 

2h"^   +  U      -   E 
m 

for 

E  > 

^n-1      ' 

P(E) 

Ci 

.V2(e^-E)-1 

for 

E  « 

% 

P'(E) 

< 

-1 

for 

all 

s 

f'(e) 

« 

-1 

for 

E 

large, 

f'(e) 

c:i 

-V^(e^-E)-2 

for 

E  ^ 

% 

Therefore,    In  each   of   the    intervals   into  which  the    points 

e        divide    the   E-axis,      F(E)      is   a   continuous   decreasing 

function  which  goes    from  positive   to    negative    values. 

Hence,    in  each    such    interval,      P(E)      must   have    one    and    only 

one   zero.      Let   these    zeros   be    denoted     E    ,      v  =   0, 1,2,  . « .,n-l, 

with     E     <e,,      e      <E     <   e    ,-,  ,      e      n<E      ,.      The    index 
o  1'        V  V  v+1*        n-1  n-1 

v      is,    therefore,    the    vibrational   quantum  niiinber.      From 
(i|.l6),    it    is    seen   that   the    correction  formula    (3*13 )   has 
the   following  properties: 


(1;.17) 


D(E) 

^=^ 

^    -    % 

D(E) 

^> 

2h"^  +  U 
m 

-   E 

D(E) 

< 

2h"^   +  U 
m 

-   E 

D(E) 

> 

2h"^   +  U 
ra 

-   E 

D(E) 

^ 

\    -E 

for 

E  ;« 

% 

for 

E 

large. 

for 

E   < 

e-L 

for 

E  > 

^n-1      ' 

for 

E   c:^ 

E, 

The  last  condition  is  a  general  property  of  the  Newton- 
Raphson  method  (see  (2.12)).   A  plot  of   F(E)   and   D(E), 
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Fig.    2.      Behavior   of     P(E)      and     D(E)      Curve; 
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derived  from  the  above  analysis,  may  be  expected  to  appear 
as  shown  in  Figure  2. 

A  serious  convergence  difficulty  is  immediately  obvious 
from  the  first  property  of   D{E).   This  indicates  not  only 

that   E  ?»-  e    leads  to  a  gross  underestimate  of  the  correc- 

V  ^ 

tlon  for  such  values  of  E,   but  that  the  smallness  of 

D(E)   cannot,  by  itself,  be  used  as  a  convergence  criteriono 

The  condition  E  ;3i- e    can  easily  be  detected  by  the 

existence  of  a  large   P  (E)   and  an  Increase  in  D(E)   from 

one  iteration  to  the  next.   Another  difficulty  which  may 

occur  is  a  jumping  from  one  branch  of  the   F(E)   curve  to 

another  on  successive  iterations.   This  can  cause  one  to 

miss  some  desired  eigenvalues  and  waste  computing  effort 

in  converging  to  eigenvalues  which  are  not  wanted  or  which 

have  already  been  computed.   Convergence  problems,  there- 

fore,  are  related  to  the  distribution  of  the  vertical 

asjmiptotes  at  e  ,   which  are  determined  by  the  crossing 

point  R  ,   and  to  the  selection  of  initial  trial  values 
m 

of  E. 

To  investigate  the  role  of  R  ,   consider  the 
^  m* 

convergence  factor  (2,13), 

ik'iQ)  Q  =  f"(e^)/2p'(e^)  , 

of  the  Newton-Raphson  method.   For  ^„  ^^   ®-„i  > 


(i4..l9)  Q  ^  (e^,  -E^) 
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V      V' 

-1 


showing  that   Q  is  large  and  convergence  is  poor  if  an 

eigenvalue   E   is  near  a  vertical  asymptote.   It  was 

pointed  out  above  that   e    is  a  value  of   E  for  which 

P   =  0   on  the  inward  or  outward  integration.   Therefore, 
in 

E     ^e    ,      is  the  situation  where,  for  the  correct  solution, 
P  ^  0.   In  other  words,  convergence  is  poor  if   R    is 
near  a  node  of  the  desired  solution.   Some  practical  means 
for  selecting  the  crossing-point,   R  ,   and  the  initial 
estimates  of   E   are  given  in  the  next  section. 


5.   Application. 

A  computer  subroutine  was  written  which,  when  given 
a  numerical  potential,  an  initial  estimate  of  E,   and 
an   e,   integrates  by  the  Numerov  method  (Equation  (2,6) ) 
and  uses  (3»15)  to  correct   E,   After   D(E)   starts 
decreasing  from  one  iteration  to  the  next,  the  convergence 
criterion  D(E)  <  e   is  applied. 

In  order  to  keep   R   as  far  as  possible  from  a  node 
of  the  solution  and  to  acquire  maximum  significance  where 
the  solution  is  large,   R   is  selected  during  the  integra- 
tion as  the  point  where  |P. |  is  largest.   Since,  for  the 
applications  considered,  this  occurs  at  the  outermost 


maximum  point  of  |P.|,  the  inward  integration  is  performed 

first  and  R   is  taken  as  the  point  at  which  |P.j  stops 
m  '  1 ' 

Increasing  with  decreasing   i.   The  description  of  the 
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behavior  of  F(E),   given  in  the  previous  section,  is 
somewhat  invalidated  since,  there,  it  was  assumed  that 
R   was  held  fixed  while  varying  E.   Instead,  the   F(E) 
given  by  the  present  program  will  behave  like   P(E)   of 
the  previous  section  in  each  interval  where  the  variation 
of  E  does  not  change   R  . 

The  data  used  here  to  demonstrate  the  method  was 
obtained  for  the  case  where   U(R)   of  equation  (2,1)  is 
a  Morse  potential  [37] 

U(R)   =  D  [1-exp  (-a(R-R  ))]^  -  D   . 

e 

Values  of  the  correct  analytic  solution  are  given  for 

comparison.   The  range   0  <  R  <  10  was  used  in  all  cases, 

A  graph  of  the  calculated  D(E)   is  given  in  Figure  3  for 

E   in  the  region  of  the  two  lowest  eigenvalues.   The 

effect  of  the  shifting  R   is  evident.   As  one  might 

expect,  the  discontinuities  are  large  in  the  region  of 

the   e  's,   where  the  solution  curves  cross  with  large 

slopes.   On  the  other  hand,  there  is  a  large  region  about 

each  eigenvalue   E   where  such  discontinuities  are 
^  V 

negligible  and  where   D(E)   is  almost  a  straight  line  of 

The  parameters  used  are:   a  =  .7112i|.8,  R  =  1.9975, 

D  =  I88.I4.355.   These  were  obtained  by  fitting  the  Morse 

potential  to  a  computed  potential  curve  for  Ht. 

Energies  are  given  in  units  of  2[i   a.u.   where   ji   is 
the  reduced  mass  of  the  nuclei. 
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slope   -1   and  is,  therefore,  a  good  estimate  of  the 
necessary  correction. 

To  show  how  the  procedure  converges  when  a  poor  first 
estimate  of  E   is  used,  the  program  was  given  a  value  of 
E  near   e, .   As  one  can  see  in  Table  2,  where  the 
successive  trial  values  are  given,  and  in  Figure  3,   D(E) 
Increases  for  a  number  of  iterations  and  then  goes  down 
rapidly.   Table  3  gives  some  of  the  eigenvalues  obtained 
for   n  =  5O5  100,  150,  and  200  points.   Exact  values  of 
the  analytic  solution  are  given  for  comparison.   Table  l\. 
contains  some  of  the  values  of  the   v  =  0   eigenf unctions 
which  were  obtained  with  50,  100,  and  200  points.   Values 
obtained  from  the  exact  solution  are  given  in  the  last 
column.   Figure  1+  contains  a  graph  of  the  potential  (5.1) 
and  the  wave  functions  obtained  for  the   v  =  0,1,  and  2 
states. 
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TABLE  2.   SUCCESSIVE  ITERATES  OBTAINED  WITH  POOR  FIRST 


ESTIMATES    OP      E^      AND      E^ 

i 

g(l) 
o 

P(E^^^) 
0 

p(l) 
^1 

F(e|^^) 

1 

-168«800  00 

-1I;70. 

-168. 500  00 

1+56.7 

2 

-169,1^50  37 

-727o9 

-166.581  98 

220.5 

3 

-170„700  20 

-350.7 

-163.292  37 

58.21 

k 

-172.90i|  1|9 

-113.0 

-160.385  38 

1,651+ 

5 

~176o702   1[^ 

-26.21; 

-160,283  89 

0OOO3 

6 

-178,621+  20 

-2,01+2 

-160.283  69 

.1550  xlO" 

-1+ 

7 

-178.797   03 

-ol786xlO"-^ 

8 

-178o798  56 

-.5185x10'^ 

9 

-178o798  57 

-ol228x  10"^ 

TABLE  3.   dependence  OF  EIGENVALUES  ON 

OP  INTEGRATION  POINTS 


n,   THE  NUMBER 


n 



E 
0 

^1 

^2 

^3 

\ 

50 

-178,81052 

-160,35850 

-11|3. 02050 

-126.83300 

-111.81091; 

100 

-178. 79921; 

-160,28781+ 

-11+2.79397 

-126,31918 

-IIO.863I1-8 

150 

-178.79866 

-160,281+28 

-11^.2.78276 

-126.?9)|)|-l 

-110.81921 

200 

-178.79857 

-160.28368 

-11+2,78090 

-126.29031 

-110.81191 

Exact 

-178.79850 

-160,28182 

-1)|2. 77990 

-126,28821+ 

-110.80832 

TABLE  1;.   DEPENDENCE  OF  WAVE  FUNCTION 


P(R),   FOR   V  =  0, 


ON     n,      THE   NUMBER   OF 

INTEGRATION 

POINTS 

R 

P(R) 

n  =  50 

n  =  100 

n  =   200 

Exact 

1.2 

,022   5903 

,022   3875 

.022  3751; 

.022   37l;6 

1.6 

.1+81+  0971; 

,1+85  88[|i; 

.1;85  9861+ 

.1;85  9927 

2.0 

1.312  2858 

1,310  2630 

1,310   11;71 

1.310  11;05 

2.1+ 

.732  9700 

«731;  7376 

.731;  81;13 

,731;  81;76 

2.8 

.126  31+1+8 

.126  1;907 

.126  1;998 

.126   50OI; 

3c2 

.008  9571+ 

,008  951;!; 

.008  951;l 

.008   95'n 
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6,   Conclusions. 

The  above  results  indicate  that  fairly  good  accuracy 
can  be  achieved  rapidly  by  the  above  procedure,  but  that 
good  first  trial  values  are  quite  important,  not  only  in 
reducing  the  number  of  iterations,  but  in  assuring  that 
no  desired  solutions  are  missed.   Good  first  estimates 
of  E    can  be  obtained  by  fitting  an  analytic  potential 
for  which  eigenvalues  are  known  and  using  it  to  obtain 
first  estimates  of  the  lowest  eigenvalues.   An  extrapola- 
tion of  calculated  eigenvalues  can  then  give  an  approxima- 
tion to  each  new  eigenvalue. 

If  it  is  inconvenient  to  obtain  first  estimates  of 
E  in  this  manner,  one  can  have  the  computer  subroutine 
perform  just  one  iteration  for  each  of  a  series  of  values 
of  E   to  get  a   D{E)   curve.   Then,  in  each  interval  where 
D(E)   changes  from  positive  to  negative,  there  is  an  eigen- 
value for  which  a  first  estimate  can  be  obtained  by 
interpolation  of  D(E).   A  count  of  the  nodes  in  the 
solution,  given  by  the  subroutine,  is  a  check  on  the  vibra- 
tional quantum  number  assigned  and  assures  that  no  desired 
solution  is  missed. 
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V      THE   EVALUATION    OF   MOLECULAR   PROPERTIES 

1.      The    Calculation   of   Energy   Levels    and   Band   Spectrum 
Constants . 
In   the   Born-Oppenheimer   approximation,    the   energy 
levels    of   the   diatomic  molecule   in  an  electronic     y~-state 
are   taken  as   the    eigenvalues   of   the  radial   Schroedinger 
equation 

(1.1)  d^S{R)/dR^      =      2[i(U(R)  -  E)    S(R) 


where 

(1.2)  U(R)      =     J(J+l)/2^LR^   +   V(R) 

(1.3)  V(R)      =     Z^Z^R"^  +  Eq-l(R)      . 

The  integer   J   is  the  rotational  quantum  number,  [i   is 
the  reduced  mass  of  the  nuclei,   ^a'^)^    ^^  ^^®  ordinary 
Coulomb  repulsion  energy  of  the  nuclei,  and  E  ,  (R)   is 
the  electronic  energy.   It  is  assumed  for  convenience  that 
V(oo)  =  0.   In  the  usual  terminology,  the  equilibrium 
distance   R  ,   the  binding  energy   D,   and  the  force 
constant   y   a-i'©  defined  in  terms  of  (1.3)  by 


[dV/dR]j^^j^   =   0 
e 

(l.i^.)  •  D  =   -  V(Rg) 

Y  =   [d^V/dR^Jj^^R 
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e 


Morse  [37 J  has  constructed  the  hypothetical  inter- 
nuclear  potential  function 

(1.5)  \iR)      =     D[l-e"^^^"^e^]2  -  D 
where 

(1.6)  a  =   -/y72D   . 

If     Vj^(R)      is    substituted  for      V(R)      in    (1.2),    the   result- 
ing energy  levels,    or   terra   values,    are   given   [38]   by 

T    ,     =      -   D  +  cJ  (v+  3)    -  /aJX    (v+  h^  +  B^J(J+1) 
vJ  e  2  e   e  2  e 

(1.7)  -      D  J^(J+1)^    -    a   (v+i)    J(J+1)      . 

e  e  d 

where 


(1.8)  u)       =      a -/2D7iI 


e 


(1.9)  Xg      =     U)Jl^ 


(1.10)  B        =      {2[i.R^)'-^ 

e  e 


(1.11)  D^      =      [43^0)^ 

e  e      e 


(1.12)  a^      =      S{Yx^bI/oJ^    -   B^/uJ^)       . 

Equations    (I.6)   and    (1.8)    give 


(1.13)  Y     =      cJq   ^i      . 


-    nil  - 


In  (1.7)  the  constant  term,  -D,  and  the  ^J  or  B^ 
terms  are  precisely  those  which  arise  in  the  case  of  a 
heirmonic  oscillator  with  force  constant  y  or  a  rigid 
rotator  with  an  internuclear  separation  R  ,   respectively. 

Energy  formulas  for  experimental  energy  levels  are 
usually  determined  by  including  higher  powers  of   (v+p) 
and   J(J+1)   and,  in  the  notation  of  the  Bohr  theory, 
they  appear  in  the  form 

(l.ll;)  T    ,      =      -   D  +   G(v)    +   F    (j) 

VJ  v 

where 

(1.15)  G(v)      =    U)   (v+^)-a;x    iv+h^  +  cJY^iv+h^   + 

+     ....  z    (v+^)^  +    ... 
e   e  2 

(1.16)  F^( J)       =      B   J(J+1)    -    D^J^(J+1)^    +   H^J^(J+1)^    +    ... 

(1.17)  B        =      B      -    a    (v+  ^)    +   Y    (v+i)^    5    (v  +  ^)^   +    ... 

V  ee2e2e2 

(1.18)  D        =     "i^     ^   ^    (v  +  ^)    +    .  ..       . 

V  e        "^e  2 

In  a  more  compact  notation,  (l.lLj.)  may  be  written 

(1.19)  T  ,  =   5~Y  .(v+i)^  jJ(j+l)J 

where  the  identification  with  the  coefficients  of  (1.15)- 
(1.18)  is  obvious. 
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Dunham  has  shown  [39],  in  his  studies  of  the  finer 
interactions  of  vibrational  and  rotational  motions,  that 
if  the  radial  internuclear  potential  (1.2)  is  of  the  form 

Uj^(R)   =   -  D  +  a^C^(l  + a-L€+  a^C^  +  ...) 

(1.20)  P 

+  B  J(J+1)  (1-  24  +  3C   -  ...  )  , 

where   4  =  (R-R  )/R  ,   the  energy  levels  will  satisfy  an 
equation  of  the  form  (1.19)  with  Y  . 's   which  he  gives  in 
terms  of   u,   and  the  coefficients  of  (1.20).   For   B  /i^ 
small,  the  band  spectrum  constants  resulting  from  Dunham's 
procedure  are  approximately  equal  to  the  ones  appearing 
in  the  Morse  solution  (1,7). 

Thus,  by  Introducing  the  effects  of  the  anharmonicity 
of  U(R),   the  energy  levels  and  band  spectrum  constants 
are  seen  to  deviate  from  those  given  by  the  ideal  relations 

D*  =   -  V(R^) 
e 


(1.21)  Jt      =   (y^)-^^^ 


e 


B*   =   (2^tRf)-l 
e         e 


satisfied  by  the  Morse  solution. 

In  Van  Vleck's  analysis  [[4.0]  of  the  dependence  of 
molecular  constants  upon   jo,   in  the  various  isotopic 
species  of  a  diatomic  molecule,  estimates  were  obtained 
of  certain  corrections  which  cause  the  ^.-dependence  of 
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experimentally-determined  values  of  D,  CJ,      and  B 

to  deviate  from  that  given  by  (1.21).   The  corrections 

considered  by  Van  Vleck  are  due  to  (a)  the  effect  of 

anharmonicity,  (b)  the  coupling  between  vibrational  and 

electronic  motion  and  the  fact  that  the  center  of  mass 

of  the  nuclei  is  not  the  same  as  the  center  of  mass  of 

the  molecule,  and  (c)  and  (d),  the  error  introduced  by  the 

separation  of  nuclear  coordinates  in  the  wave  function  of 

the  molecule.   He  derived  expressions  for  the  (b),  (c) 

and  (d)  corrections  to  cJ   and  B    and  estimated  them 

e        e 

for  the  case  of   Hp,   HD,   and   Dp  by  using  the  Wang 
wave  functi,on.   He  derived  a  theoretical  estimate  of  the 
error  resulting  from  the  assumption  of  infinite  nuclear 
masses  which,  with  known  corrections  for  the  separated 
atoms,  can  be  used  to  correct   D.   He  expressed  the 
corrections  to  the  "ideal"  quantities  (1.21)  in  the  form 


D   =   d""  +  A.  +  Ac 


(1.22)  ^^        =    uJ^  (1+  5„+  5 J 

e      e      as 

B    =   B^*  ( 1  +  Y  +  Y  )   • 
e      e      '  a   '  s 

The  subscript   "a"   refers  to  the  anharmonicity  correction 
and  the  subscript  "s"   refers  to  the  sum  of  the  (b)-(d) 
corrections.   In  the  present  work,  the  term  values  are 
found  by  solving  (l.l)  by  the  ni:imerical  method  of 
Chapter  IV.   Consequently,  the  molecular  constants  derived 
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from  them,  which  will  be  distinguished  by  the  super- 
script "-K-!!-"  ,   need  only  be  corrected  for  the  (b)-(d) 
effects.   In  Van  Vleck's  treatment,  the  correction 
formulas  may  be  written 

D  =   D""  +  A3 

(1.23)         a)g  =  oJ^   (1+  5g) 

e      e       s 

For  the  purpose  of  comparing  the  various  corrections  to 
the  ideal  quantities  (1.21),  an  estimate  of  the  anharraoni- 
city  correction  may  be  derived  from  the  relations 

(1.21;)  .j'"'  =  4'  (1+5J 

B*-"^   =  b''  (1  +  Y  )   . 
©  ■       e      a 

Before  forming  quantitative  relations  between  experi- 
mental data  and  an  internuclear  potential,  some  considerations 
of  the  magnitude  of  the  uncertainty  of  each  quantity  involved 
must  be  made.   These  are: 

(1)  Van  Vleck's  corrections  (b),  (c),  and  (d),  are  based 
on  estimates  of  their  effects  on  the  harmonic  oscillator  and 
rigid  rotator  and,  consequently,  cannot  be  applied  to  quanti- 
ties resulting  from  motions  departing  too  greatly  from  these 
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idealized  situations. 

(2)  The  values  of   D  ,   B  ,   and  particularly  '^ 

e  e 

are  quite  sensitive  to  both  random  and  systematic  errors 
in   V(R). 

(3)  Calculated  term  values   T  j      are  affected  by  the 
systematic  errors  in   V(R)   which  result  from  the  fact 
that  the  accuracy  of  trial  electronic  wave  functions  may 
vary  over  the  range  of   R.   F\arthermore ,  with  a  fixed 
number  of  values  of   V{R),   the  numerical  accuracy  of 

T' '   decreases  with  increasing  v. 

(il)  There  is  a  certain  amount  of  ambiguity  in  the 
definition  of  the  band  spectrum  constants.   In  practice 
one  finas  that  the  polynomials  in  v  + -p   and  J(J+1) 
which  one  tries  to  fit  to  the  vibrational  and  rotational 
levels  alternate  in  sign  and,  consequently,  their  coeffi- 
cients depend,  quite  critically,  upon  the  nxomber  of  term 
values   T  y   used,  the  type  of  litting  (i.e.,  exact,  least 
squares,  etc.),  and  the  accuracy  of  the  term  values.   An 
even  more  basic  difficulty  is  that  low  order  polynomials 
in   v+  ^   and   J(J+1)   having  coefficients  ^^        and  B 
which  can  be  used  to  determine   y   and   R    simply  fail 
to  represent  the  totality  of  term  values. 

A  brief  perusal  of  the  derivation  of  band  spectrum 
constants  from  experimental  and  theoretical  data  tends  to 
indicate  that  in  the  calculation  of  the  coefficients  of 
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F  (J),   (1,16),  most  of  the  indeterminacy  is  caused  by 

the  rapid  loss  of  significant  digits  in  the  successive 

coefficients  rather  than  a  failure  of  the  energy  formula 

(1.16).   Results  also  indicate  that  the  effect  of  the 

ambiguity  in  B   is  small  relative  to  both  the  Van  Vleck 

correction  to  B   and  the  numerical  and  experimental 

e 

errors.   The  difficulties  arise  when  one  tries  to  fit 
polynomials  in  v+ -^  to  the  values  of   G(v),  B  ,  D  ,  o.o 
or  their  differences.   An  unambiguous  comparison  with 
experiment  can  be  made,  however,  by  considering  the 
vibrational  quanta 

(1.25)  V  G(v)   =   G(v)  -  G(v-l)   =  6j  -  2&y  X  V  +  o..  . 

e    e  e 

For  the  lowest  vibrational  states,  to        is  much  larger 
than  the  sum  of  all  other  terms  in  (1.25)  so  that  one 
can,  within  the  framework  of  Van  Vleck 's  assumptions, 
apply  his  correction  to  the  vibrational  quanta  S/   G(v) 
instead  of  cj       and  obtain  a  more  direct  comparison  with 
observed  spectra.   Fiorthermore,  since  the  correction  is 
in  the  form  of  a  multiplicative  factor,  a  properly 
corrected  60   will  result  no  matter  what  type  of  curve- 
fitting  is  used  to  determine  the  coefficients  of  (1.25) 
from  the  corrected  vibrational  quanta.   By  the  same 
argument,  one  may  also  apply  Van  Vleck's   B    correction 
factor  to  the  rotational  term  values   F  (j)   and  to  the 
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rotational  constants   B  .   These  considerations  form  the  basis 
for  the  procedure  used  in  Chapter  VII  for  comparing  the 
calculated  and  experimental  spectral  data  for  Hp,  HD,  and  Dpo 

2.   Averaging:  over  R« 

After  having  obtained  numerical  solutions   S(R),  T(R), 
...   of  the  radial  Schroedinger  equation,  one  may  proceed  to 
calculate  expectation  values  and  transition  matrix  elements 

00 

(2.1)         Q  =  J    S(R)  Q(R)  T(R)  dR 

b 

corresponding  to  an  operator   Q(R).   Denoting  the  integrand 

by   f(R)   and  assuming  that  ^^r+l      ^^    ^°  large  that   f(R) 

and  all  of  its  derivatives  are  negligible  at   R  =  R^+n > 

the  -ciulerian  sum   formula   gives 

(2.2)  Q  .,rf(R)dR  =  hg:f,+|f^-|^f; +^f';-.. 


where 


^1   ^  ^^^i^  • 


In  most  cases  the  solutions   S(R),  T(R),  ...   and  the 
operator   Q(R)   are  such  that   f(R)   and  its  derivatives  to 
high  order  are  zero  at   R  =  0.   Therefore,  the  trapezoidal 
rule,  which  results  from  neglecting  derivatives  in  (2.2), 
is  an  integration  formula  of  quite  high  order  accuracy 
in  h. 

"■  Equation  (5.8.20),  Page  lSl\.,    Ref.  [I4.IJ. 
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VI   AN  LCAO  CALCULATION  FOR  E^ 

1.   Introduction. 

It  is  immediately  obvious,  in  a  calculation  of  this 
type,  that  it  would  be  desirable  to  vary  as  many  wave 
function  parameters  as  possible  to  minimize  the  electronic 
energy.   However,  for  many-electron  molecules,  the  comput- 
ing effort  involved  in  obtaining  all  the  integrals 
required  in  a  full  variation  of  the  orbital  exponents  of 
the   AO's  would  make  the  cost  of  such  computations 
prohibitive.   On  the  other  hand,  as  one  can  see  from  the 
form  of  the   STO's   (II  L|..  1 ) ,  applying  a  scale  factor  y? 
to  all  orbital  exponents  is  equivalent  to  applying  the 
scale  factor  to  all  electronic  distances  in  the  molecule. 
The  special  role  of  the  parameter  >7  ^^"^  methods  for 
optimizing  it  have  been  discussed  by  LBwdin  [I4.2J.   Such 
methods  are  particularly  convenient  since  expectation 
values  of  the  potential  and  kinetic  energies  of  the 
electrons  and  the  Integrals  used  to  compute  them  satisfy 

(1.1)  V(R,  yi  )      =      ^  V(R  yi,l) 

and 

2 

(1.2)  T(R,  >^  )  =     \    T{R  y^,!)      , 

respectively.   This  means,  in  effect,  that  a  table  giving 
the  integrals  as  functions  of   R,   with   ^j  =  1,   can  be 
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used   to    obtain  the    integrals    for   arbitrary      R      and      >;;   <, 
Furthermore,    the   same   tabular   values   used  to   interpolate 
the    integrals   can  also  be   used   to  evaluate   the   R-derivatives 

of      V      and      T      in  LBwdin's    formula   for      ^)E/d,;. 

t 

The   Hp  molecule,  which  does  not  require  the  compli- 
cated two-electron  integrals,  is  an  excellent  test  case 
for  studying  the  role  of  the  orbital  exponents  and  for 
comparing  the  results  of  varying  a  scale  factor  with  the 
results  of  varying  the  orbital  exponents  independently. 

In  addition  to  the  above  reasons  for  carrying  out  a 
calculation  of  approximate  solutions  where  exact  analytic 
solutions  have  already  been  obtained  [I4.3],  there  is  some 
interest  in  having  theoretical  values  for  certain  quanti- 
ties which  can  be  calculated  with  these  approximate 
solutions  and  the  machine  programs  written  during  the 
course  of  this  work. 

2.   Solution  of  the  Electronic  Wave  Equation. 

In  the  present  calculations,  the  electronic  ground 

2x + 

state   2 ^^s  considered  and  the  Is,  2s,  and  2pcr 

O 

AG's   on  the  two  centers  were  used.   Due  to  the  symmetry 

of  the  molecule,  the  resoective  orbital  exponents 

u-T  ,  M-o  »   and   ii,„    were  the  same  on  each  center  and 
^Is'  ^2s'        ^2pcr 

there  were  only  three  g-type  MO's  to  consider.   The  MO's 
were  calculated  by  the  method  described  in  Section  II. Ij. 
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and   the      3x3      matrix      H      of   the   Hainlltonian 

(2.1)  'M  =   -\/\-^-^  -r-^ 

with  respect  to  the  g-type  MO's  was  diagonalized. 

The  calculation  was  programmed  as  a  subroutine  which, 
for  a  given  internuclear  distance   R   and  a  given  set  of 
three  orbital  exponents,   ix-,  ,  [Xp  ,   and   jXp   .  ,   yielded 
three   ^      wave  functions  and  energies.   For  the  rest 

of  the  discussion,  the  lowest  of  these  energies,  which  is 

+ 
an  approximation  to  the  electronic  ground  state  of   Hp, 

will  be  referred  to  as  a  function  ^  ^^jI^io' t^Ps' '^2Dcr  ^   °"^ 
the  four  parameters  given  to  the  subroutine. 


3.   Internuclear  Force  Curve. 

For  an  exact  solution  of  the  electronic  wave  equation, 
the  internuclear  force  due  to  the  electrons  in  a  molecule 
is  given  by  the  Virial  Theorem  (VT) 

(3.1)  dE/dR  =   -  (2T+  V)/R 

where   T,  V,  and  E   are  the  kinetic,  potential,  and  total 
energies,  respectively,  of  the  electrons.   Under  the  same 
conditions,  the  internuclear  force  is  also  given  by  the 
Hellman-Peynman  Theorem  (HFT), 


/  N    p        \ 
(3.2)      dE/dR   =   (  |II  r^^  cos  9^^^ 


\i=l  °-^  "y  av.el. 
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The  right  side  of  (3.2)  represents  the  average  value,  over 
electronic  coordinates,  of  the  electric  field  at  nucleus 
A   due  to  the  electronic  charge  cloud.   It  can  easily  be 
shown*  that  the  VT  will  be  satisfied  by  a  wave  function 
for  which  a  scale  factor,  applied  to  electronic  coordinates, 
minimizes   E   at  each  R.   Here,  and  in  what  follows, 
"minimize  E"  will  imply  that,  at  the  value  of  the  parameter 
in  question,   E   is  a  continuous  dif f erentiable  function 
of  the  parameter  and  has  a  relative,  but  not  necessarily 
an  absolute,  minimum.   Hurley  [L|i4.]  has  discussed  the  condi- 
tions under  which  approximate  wave  functions  satisfy  both 
the  VT  and  the  HFT.   Hurley's  results,  in  the  present 
application,  show  that  the  HFT  will  be  satisfied  if  all 
parameters  in  the  electronic  wave  function  which  are 
allowed  to  vary  with  R  minimize   E   at  each   R.   The 
linear  coefficients  of  the  MO's  and  optimized  orbital 
exponents  satisfy  these  requirements.   The  one  remaining 
R-dependent  parameter  in  the  present  wave  function  gives 
the  locations  of  the  centers  of  the  AG's.   Hurley's 
results  imply  that  the  HFT  would  have  been  satisfied  if 
the  calculation  had  been  performed  with  floating  orbitals, 
i.e.,  orbitals  having  centers,  not  at  the  nuclei,  but  at 
symmetric  points  on  the  nuclear  axis  determined  by  minimiz- 
ing E  with  respect  to  the  distance  between  them.   (Use  of 

''^  See  LBwdin,  Ref.  [i|2j. 
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floating  orbltals  would  necessitate  the  use  of  some  three- 
center  nuclear  attraction  integrals  and  would  have  intro- 
duced one  more  nonlinear  parameter  in  the  variation  process.) 

The  variation  of  orbital  exponents  was  performed  in 
two  ways.   In  the  first,  which  will  be  called  Case  A,  all 
the  orbital  exponents  were  set  equal  to   >^  .   At  each  of 
the  internuclear  distances   R  =  ,125(. 125)10  one  or  two 
values  of  r]      which  minimize   E{R,  )^  ,  ^  ,  >7  )   were  found. 
Letting  Y[  (R)   denote  the  value  of  V)      which  gives  the 
lower  minimum  of  E,   >o  (R)   was  found  to  jump  from  one 
continuous  curve  to  another  between  R  =  I4..250   and 
R  =  I4..375.   The  Y^  (R)   curve  Is  plotted  in  Figure  5  with 
nonconnected  points  indicating  values  of  h   where  the 
higher  minima  were  found. 

The  results  of  the  calculations  indicate  that  if  one 
optimizes  )|   at  a  relatively  few  points  and  Interpolates 
to  get  values  of  vi    (R)   at  other  points,  one  must  optimize 
^l      at  enough  values  of  R   to  establish  that  the  points 
used  In  each  application  of  the  Interpolation  formula  lie 
on  the  same  continuous  curve.   One  further  consequence  of 
the  discontinuous  )^  (R)   curve  is  that  the  VT  and  HFT  will 
not  hold  at  points  of  discontinuity  In   m  (R)   or  at  points 
obtained  by  Interpolating  across  such  discontinuities. 

In  the  second  calculation,  which  will  be  called  Case  B, 
the  three  orbtial  exponents   M-Toj  ^2s'      ^^^      '^2Dcr  ^^^^ 
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Pig.  5,   Graph  of  Optimal  Scale  Factor  ,yj{R) 

varied  Independently  by  the  procedixre  described  in  Section  II 
10.   This  added  the  following  complications  to  the  problem: 

1.  There  were  more  minima  at  each  R.   For  example,  six 
minima  were  found  at  R  =  2.0. 

2.  More  points  were  found  where  the  best  minimum  changed 
from  one  continuous  curve  in  the  (R,ti-,  ,[i^    ,[i^        ) -space 

to  another. 

3.  The  energy  surface  defining   E   as  a  function  of  the 

M-'s   for  fixed   R  was  often  quite  flat  in  a  diagonal 

direction  in  the  M'ig"M-2s   plane.   For  this  reason,  fewer 

significant  digits  could  be  calculated  for   a,   and  u.^  . 

Is      '^2s 
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The  Implication  of  these  results  is  that  so  little 
significance  can  be  attached  to  interpolated  values  of 
optimal   |j.'s   that  it  is  of  dubious  value  to  spend  much 
computing  effort  on  optimizing  the   ^x '  s   at  a  few  points 
and  then  interpolating  to  obtain  the  rest.   Instead,  it 
would  probably  be  better  to  optimize  the   ^.'s   independ- 
ently at  one  point  near  R  ,   the  equilibrium  value  of  R, 
and  optimize  only  a  scale  factor  on  the  resulting   ^l's   at 
all  other  values  of  R, 

In  order  to  check  the  numerical  results  and  to  make 
one  further  point,  the  numerical  derivative, 

E^   =   (E(R+.05)  -  E(R  -  .05))/.l   , 

of  the   E(R)   of  Case  B  was  calculated  and  compared  with 
the  values  given  by  formulas  (3.1)  and  (3o2).   The  latter 
will  be  denoted  by  E^^p   and  E^^m,   respectively.   Some 
of  the  results  are  given  in  Table  5«   It  was  found  that  the 
discrepancy  between  E^   and  E^m  was  less  than  the  error 
of  the  niXTfierical  differentiation  formula.   The  numerical 
results  obtained  here  also  provided  some  interesting 
information  about  the  validity  of  the  assumption,  made  by 
some  workers,  that  the  force  constant  can  be  evaluated  by 
supposing  that  the  electronic  charge  cloud  is  approximately 
stationary  for   R   near  the  equilibrium  distance   R  . 
Under  this  assumption,  one  can  differentiate  the  expectation 
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TABLE  5.   ELECTRONIC  ENERGY  OF   H    AND  ITS  DERIVATIVES 


1     R 

1 

E^^(R) 

^N 

^1 

/T 

t 

^HFT 

^V^R^ 

%.^^^ 

i    .25 

-1.398  i;85 

.619 

602 

.622 

391 

.621 

515 

.751    61; 

.881    1;12 

.50 

-1.231;  853 

.61;8 

126 

.61;8 

659 

o61;7 

820 

-.239    23 

.751;  080 

.75 

-1.002   259 

.566 

998 

.567 

li>0 

.567 

322 

-.371    39 

.571    2hh 

1.00 

-.951  662 

.1^79 

176 

.1^79 

169 

.1;79 

91;5 

-.333  75 

.1;28   079 

1.25 

-.81^1  639 

.14-03 

I4.05 

.1^03 

215 

.1+05 

621; 

-.276   02 

.325  528 

1.50 

-.714-8  765 

.31|1 

830 

.31|1 

777 

o31;3 

183 

-.220   01 

.21;8  91;6 

1.75 

-.66V   813 

.291 

661 

.291 

588 

.293 

210 

-.180   [;0 

.1V3   306 

2,00 

-.602  199 

.250 

70I4. 

.250 

759 

.252 

522 

-.11+7  58 

.152   152 

2,ZS 

-.51^3  871 

.217 

068 

.217 

026 

.219 

010; 

-.121   61; 

.121   281; 

2.50 

-.1^93  223 

.189 

767 

.189 

01;2 

.191 

191 

-.102   39 

0O97  586 

2.75 

-.lil|8  903 

.165  SSS 

.165 

583 

.167 

662 

-.086   70 

.078   820 

3.00 

-.J|"IO  072 

.11^5 

871 

.11;5  559 

.11;7 

622 

-.071;  93 

.061;  232 

3.25 

-.375  883 

.128 

361; 

.128 

31+9 

.130 

737 

-.061;  71; 

.052  751; 

3.50 

-.314-5  685 

.113 

5i;7 

.113 

529 

.116 

029 

-.055  06 

.01;3  1;91; 

i|.00 

-.295  236 

.089 

337 

.089 

361; 

.091 

871 

-.ol;2  31; 

.029   911; 

5.00 

-.223   829 

.056  31|.0 

.056 

ISS 

.058 

722 

-.025  33 

.011;  762 
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value  of  the  electronic  energy  twice,  with  nucleus  B  . 
held  fixed,  to  get 

(3.3) 

Values  of  c3  E/^R  ,  obtained  from  central  differences  of  E^rp, 
and  Q   (R),  the  electric  field  gradient  along  the  molecular 
axis  at  nucleus  A,  are  listed  in  Table  5  where  it  may  be 
observed  that  c)^E/t)R^   and   -Q^^(R),   although  quite  differ- 
ent  for  most  values  of  R,   differ  by  only  3 /o  at  R  =  R  =  2.00. 
It  has  been  ghown  by  W.  Clinton  \.hS\    that  the  assumption  of 
a  static  charge  cloud,  taken  with  the  VT  and  HPT,  leads  to 
contradictions  and,  is  therefore  invalid.   He  also  points 
out  that  deceptively  good  agreement,  based  on  this  false 
assumption,  does  seem  to  occur  with  the  use  of  LCAO  wave 
functions.   This  observation  is  supported  by  the  present 
calculation  where  there  is  a  discrepancy  of  only  3  h  between 
-  d^^/d^      and  Q   (R)   at   R  =  2.00  while,  in  the  more 
accurate  non-LCAO  calculation  of  Auffray  and  Gooley  \.\\h] 
Q^^(2)   is  18%  greater  than   -  [^^^E/^R^ Jj^^2  " 

The  dissociation  energy  in  Case  A  was  .102  185  a.Uo, 
for  Case  B  it  was  .102  199  a.u.,  and  in  the  accurate  calcxiLa- 
tion  of  Auffray  and  Gooley,  it  was  .102  63i|-  a.u.   For  other 
values  of  R,   the  difference  between  Case  A  and  Case  B  was 
slightly  greater  but  hardly  enough  to  warrant  the  amount  of 
calculation  involved  in  varying  three  nonlinear  parameters 

instead  of  one, 
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VII   CALCULATIONS  FOR  E^,      HD,   AND   D^ 

1,   The  Electronic  Wave  Functions  and  Internuclear  Potentials. 

The  results  reported  In  this  chapter  were  calculated 

with  two  different  electronic  wave  functions  for  the   > 

g 

state.      The   first   of   these  was    calculated  by  McLean,    Weiss, 
and  Yoshlmlne    [l\.7 ]    (hereafter  referred  to  as  MWY)    and  was 
of   the   form 

(1.1)  1     =     H   C      I 

k=l     ^     ^ 

where 

ll     =      ^^^g'^^g^    '      ^2     "      ^2Sg,2p<r  )    , 

(1.2)  I3      =      (13^,13^)    ,      \     =      (2p7i^,2p7t^), 

i^     =      (2pitg,2p7rg)    . 

Since  an  antisymmetric  spin  function  can  conveniently  be 
factored  out,  ^  Is  constructed  as  a  symmetric  function 
of  space  coordinates.    The  notation  In  (1,2)  denotes  the 

symmetrized  products 

1  I  I 

(1.3)  i^,X    )      =     X(l)     /  (2)  +  X{2)     /(I) 

where  the  one-electron  space  functions  are  nonorthonormalized 

primitive  symmetry  orbltals, 

_1 

(l.i;)  X  =   2*2  (/  +/  )  , 

^  '^  g  or  u         ^a  —  b   * 
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constructed  of  the  AO's  In  (II  I]..l).   The  STO's   Is   and 
Is    in  J,   differ  only  in  the  orbital  exponent.   The  five 

orbital  exponents  [i.^^^,    l^^^g,,  t^js'  ^2po- »  ^^^     ^2p7t  ^®^® 
optimized  at  five  internuclear  distances  and  their  values 
at  the  remaining  points  were  obtained  by  interpolation. 
WdY   state  that  a  minimum  was  assumed  to  have  been  found 
when  three  significant  figures  agreed  in  the  binding 
energy  on  successive  iterations.   However,  in  the  minimiz- 
ing for   Hp  described  in  the  previous  chapter,  where  seven 
significant  figures  was  the  criterion,  there  were  six 
minima,  at   R  =  R  ,   all  agreeing  in  the  first  three  digits 
but  having  quite  different  orbital  exponents.   Therefore, 
MWY's  orbital  exponents  are  possibly  not  the  ones  for  the 
absolute  minima.  This,  and  the  possibility  of  a  discontinu- 
ous optimal  orbital  exponent  curve  lead  one  to  doubt  the 
accuracy  of  the  interpolated  orbital  exponents.   The 
present  author  would  suggest  the  possibility  that  a  more 
thorough  calculation  of  optimal  orbital  exponents  near  the 
equilibrium  distance,  followed  by  a  variation  of  a  scale 
factor  for  other  values  of   R  would  have  given  a  better 
internuclear  potential  curve  with  about  the  same  amount  of 
effort. 

With  the  orbital  exponents  so  determined,  MWY  calculated 
the  coefficients   C,   of  (1.1)  at  1^5  non-unif ormly  spaced 
points  in  the  range   .8  <  R  <  6,0. 
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In  the   present    calculation,    the    values   of 

1~1 


(1.5)      W^^^R)      =    /l*  IZ  r'J"^   "'    fn^°°^   ®al^   ^  "^^   ' 


where   P    is  the  n-th  Legendre  polynomial,  were  obtained 
for   n  =  0,  1  and  2  at  each  value  of   R  for  which  MWY 
gave  the  wave  function  parameters.   The  program  which 
calculated  these  average  values  was  used  to  calculate  the 
normalization  factor  as  a  check.   Then,  the  values  of 
(1.5)  and  the  electronic  energy  E  ,(R)   were  obtained  at 
the  equally-spaced  points   R  =  0(.02)6   by  quadratic 
interpolation.   This  permitted  the  use  of  the  program 
described  in  Chapter  IV  to  calculate  the  nuclear  wave 
functions.   The  spacing   h  =  .02   is  that  given  by  MWY 
around  the  equilibriiJm  value  of   R.   The  maximum  R  given 
by  MWY  was  found,  by  numerical  trials,  to  be  high  enough 
to  solve  the  radial  Schroedinger  equation.   For  the  low 
range  of   R,   the  two  lowest  calculated  points,   R  =  .8 
and   R  =  .9,   and  values  at   R  =  0   appropriate  to  the 
united  atom,   H  ,   were  used  by  the  interpolation  proced- 
ure.  Errors  introduced  in  the  internuclear  potential  for 
low  R  were  relatively  small  due  to  the  large  Coulomb 
repulsion  term  1/R   and  would  be  of  importance  only  if  one 
were  to  try  to  go  to  higher  vibrational  levels. 

The  second  electronic  wave  function  was  calculated  by 
Kolos  and  Roothaan  [L|.8J  (hereafter  referred  to  as  KR)  and 
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was   of  the   form 

(1.6)  1     =    IZ  ^i\ 

where 

|(p,q,r,s,^.)      =      6  ^2      ^p^q   ^r^s   ^^-^ 

Here,      £,     and     )^      are   the   elliptic   coordinates   and      a     is 
a   scale   factor  which  was   optimized  at   each  of  the    inter- 
nuclear   distances   at   which   the    calculation  was   performed. 
Except   for    large      R,      the   KR   electronic    energies   are 
lower   than  those   obtained  with   the   MWY  wave   function.      The 
internuclear   potential    for   the    full   range    of      R,      at 
intervals   of      .05   a.u.,      was   obtained  by  using  MWY's   values 
for   high     R      and  by  interpolating  as   described  above.    The 
equilibrium  energies    for   the  MWY   and  KR  potentials  were 
36   6l\hr.09   cm."-*-      and      38   285.60   cm."-*-,      respectively. 

2.      Vibrational-Rotational  SpectriJm, 

A   large   number   of  vibrational-rotational  wave   functions 
and   term  values   for     Hp,    HD,    and  D^     were    calculated  with 
both  the  MWY   and   the   KR   internuclear   potentials   using   the 
reduced  masses 

^l     =     918.06  IJ-'      =     1223.77  [^       =     l831^.709l|. 
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for      Hp,    HD,    and   Dp,      respectively.       Some    of    the    term 
values    obtained  with   the   more    accurate   KR   potential    are 
listed   in  Table    6  with    the    zero-point    taken   as    the    energy 
( -1    a.u. )    of   the    dissociated   atoms,      H(ls)    +   H(ls). 

The    correctiohs    estimated   by   Van   Vleck  for      Hp,    HD, 
and  Dp,      which  will   be   used  here,    are,    in   the   notation  of 
(V  1.23), 

As      =     ^'5   cm"^        A's      =     h^e   cm'l       Ag   =   2 . 8   cm"l 


(2.1)       5         =      -.00071  5'       =      -.00052  6^    =    -.00035 


I 

's  ' "  "s  ""  's 


II 


Y         =      -.00089  y'       =      -.00068  Y      =    -.000i;i4. 

s  s  s 

where  the  unprimed,  primed,  and  double-primed  quantities 
are  for   Hp,  HD,  and  D^,      respectively. 

The  vibrational  quanta  (1,25)  obtained  by  usinp-  the 
jWY  and  KR  potentials  are  listed  in  Table  ?.   The  result 
of  applying  the  Van  Vleck  corrections  to  the  vibrational 
quanta  obtained  with  the  KR  potential  are  given  in  column 
four  for  comparison  with  the  experimental  values  in 
column  five.   It  may  be  observed  that  the   ^fWY  vibrational 
quanta  are  too  small  by  an  amount  which  increases  almost 
linearly  in  magnitude  from  12cm    at   v  =  1   to  53  cm 
at   V  =  6  while  the  KR  results  are  too  large  by  an  amount 
xirhich  increases  only  slightly  from   8  cm     at   v  -  1. 
For     V  >  6,   the  discrepancy  in  botn  cases  is  in  the  same 
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TABLE  6.       CALCULATED      TERM   VALUES      ( -T    j) ,      BELOW  DISSOCIA- 
TION   LIMIT,    FOR      ti^,    HD,    and   D2    (IN   CM""^). 


V     ° 

1 

2 

3 

0 

36  103.50 

35   981;.  89 

35  71^8.71; 

35  397.21 

1 

31  93I1.21 

31  821. 1;9 

31  597.10 

31  263.11 

2 

28  001.11; 

27  891;.  18 

27  681.26 

j 

27  361;.  1;0 

H2  3 

21;  296.8[^. 

21;  195.55 

23  993.93 

23  693.93 

k 

20  816.97 

20  721.21 

20  530. 61; 

20  21;7.16 

5 

17  561.38 

17  1;71.18 

17  291.69 

17  021;.  73 

6 

11;  527.22 

11;  l|)|2.61; 

11;  27l;.39 

11;  021;. 22 

0 

36  393.15 

36  303.63 

36  125.80 

35  860.28 

1 

32  751;.  17 

32  668.69 

32  1;98.33 

32  210;.  27 

2 

29  293.92 

29  212.20 

29  01;9.35 

28  806.51 

HD  3 

26  006.28 

25  928.29 

2S   772.85 

25   51;l.l0 

k 

22  889.35 

22  811;. 95 

22  666.68 

22  lU^5.6i; 

5 

19  9I|1  .07 

19  870.32 

19  729.31; 

19  519.19 

6 

17  158.11 

17  090.95 

16  957.18 

16  757.82 

0 

36  737.72 

36  677.87 

36  558. 1;7 

36  380.01; 

1 

33  738.70 

33  680.95 

33  565.71; 

33   393.58 

2 

30  860.53 

30  801;.  86 

30  693.78 

30  527.82 

D2  3 

28  097.83 

28  0||)|.20 

27  937.19 

27  777.31 

i^ 

25  1;53.09 

2S  koi.k3 

25  298.38 

25  11;!;.  lA 

5 

22  918.88 

22  869.25 

22  770.22 

22   622.29 

6 

20  1;99.70 

20  1;51.95 

20  356.68 

20  211;.38 

•^Calculated  with  Kolos   and  Roothaan's   internuclear  potential 
[1;8J. 
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TABLE  7.   CALCULATED  AND  OBSERVED  VIBRATIONAL  QUANTA 
V  Giv)      FOR   H^,  HD,  AND  D^   (IN  CM"-*-). 


MWY"" 

KR^ 

KR,    Corr. 

Exp . 

1 

1+11|9.20 

1^169.29 

1 
1^166.33 

1+161.13° 

2 

3908.11 

3933.07 

3930.28 

3925.95       ' 

3 

3670,38 

37OI+.30 

3701.67 

3695.21+      1 

3l+3i|.2i^. 

31+79.87 

31+77.1+0 

31+68.01       1 

5 

3198.55 

3255.59 

3253.28 

321^1.56 

6 

2959.614. 

303I+.I6 

3032.01 

3013.73 

1 

3622.114. 

3638,98 

3637.09 

3632. 11+"^ 

2 

314^-0.77 

31+60.25 

31+58.1+5 

31+51+.  73 

3 

HD 

1+ 

3261,20 
30814.. I4.6 

3237.61+ 
3116.93 

3285.93 
3115.31 

3280.68 
3109.31 

5 

2906.80 

291^8.28 

291+6.75 

6 

2730.17 

2782.96 

2781.51 

1 

2985.71+ 

2999.02 

2997.97 

2993.51+8® 

2 

286I4..53 

2878.17 

2877cl6 

3 

271+3.36 

2762.70 

2761,73 

2625.57 

261+1+.71+ 

261+3.81 

5 

2506.I4.8 

253'!+.  21 

2S33.32 

6 

21+19.18 

21+18,33 

^Calculated  with  McLean,  Weiss,  and  Yoshimine ' s  internuclear 

potential  [1+7  J. 
^Calculated  with  Kolos  and  Roothaan's  internuclear 

potential  [1+8] . 
cHerzberg  and  Howe  [51  J* 
^Durie  and  Herzberg  [50 J. 
estoicheff  [1+9J. 
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direction  and  increases  in  magnitude  with  v.   This  can 
partly  be  attributed  to  the  nature  of  the  electronic  wave 
functions  used  by  MWY  and  KR.   The  LCAO  wave  functions  of 
MWY,  for  obvious  reasons,  give  better  electronic  energies 
at  the  extreme  values  of   R  but  yield  an  internuclear 
potential  which  is  too  shallow  in  the  intermediate  range. 
On  the  other  hand,  the  KR  electronic  wave  function  is 
better  for  the  intermediate  range  of   R   and,  consequently, 
curves  upward  too  much  for   R  very  different  from  R  . 
The  effect  is  as  one  may  expect:   the  low  curvature  of  the 
MWY  potential  causes  the  vibrational  levels  to  lie  too 
close  together  and  the  high  curvature  of  the  KR  potential 
causes  them  to  lie  too  far  apart.   However,  one  would  not 
expect  the  error  in  curvature  of  the  KR  potential  to  be 
significant  for  the  lowest  vibrational  states. 

An  unambiguous  test  of  the  theoretically-calculated 
rotational  effects  and  of  the  value  of   R    given  by  the 
KR  potential  can  be  made  by  a  comparison  with  the  pure 
rotation  Raman  S-lines  measiired  by  Stoicheff  [[4.9].   These 
are  given  by  the  differences 

(2.1)         S^(J)   =  F  (J+  2)  -  F^(J) 

In  Table  8,  Stcicheff's  values,  which  are  reportedly 
acc\irate  to   +  .02  cm.   ,  are  listed  with  the  results  of 
applying  the  Van  Vleck  correction  to  the  differences  (2.1) 
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obtained  from  the  KR  potential. 


TABLE  8.   CALCULATED^  AND  OBSERVED^  RAMAN   Sq(J)   LII^S 
OF  H^,  HD,  AND  D^. 


] 

h 

] 

iD 

^2 

Calc» 

Exp. 

% 

Dif 

Calc. 

Exp, 

% 

Dif 

Calc, 

Exp. 

7o 

Dif 

S„(0) 

35i|.l^i; 

35)+.  38 

,02 

267.17 

267.09 

.03 

179.17 

179.06 

.06 

S„(l) 

587.16 

587.06 

,02 

iUi.3.25 

i|i+3.08 

Dk 

297.70 

297.52 

.06 

S„(2) 

811;, 59 

3l[;.i;l 

J32 

616,33 

616.09 

,0i| 

l+ll4.,89 

1+11].,  66 

j06 

S„(3) 

1031^.91 

1031^,65 

^3 

785.32 

78I4..99 

JJk 

530.21 

52V. 91 

M 

s„{U) 

121^6. k2 

91^9  c  19 

9i;8.82 

.Ok 

61^3.20 

61j.2.8l 

/)6 

s„(5) 

Ikkl . 69 

1107.06 

753.39 

Calculated  with  KR  potential, 
has  been  applied. 

'Stoicheff  [l;9Jo 


The  Van  Vleck  correction 


The  rotational  constants   B  ,  D  ,  and  H  ,   eiven  in 
Tables  9  and  10,  were  calculated  by  an  exact  fit  of  the 
rotational  energy  formula  (V  1.16)  to  the   J  =  0,1,2,  and  3 
rotational  levels.   It  was  found  that  using  one  or  more  or 
less  J-levels  affected  B   only  in  the  last  figure  given. 
This  and  the  fact  that  the  percentage  deviations  from 
experiment  in  the  Raman  S-lines  of  Table  8  are  quite  close 
to  those  for  the   B^  values  of  Table  9  indicate  that  the 
uncertainty  due  to  the  procedures  for  obtaining  the   B  's 
is  negligibly  small. 
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TABLE  9.   CALCULATED  AND  OBSERVED  ROTATIONAL  CONSTANTS 


B^   FOR  E^,    HD, 


AND   D2   (IN  CM"-^). 


V 

M¥Y^ 

KR^ 

KR,    Corr. 

Exp. 

0 

59.283 

59.398 

59.31^5 

59.339° 

1 

56,271 

56.1+1+9 

56.399 

56.369 

2 

53.310 

53.567 

53.519 

53.k7$ 

H^        3 

50.390 

50.732 

50.687 

50.626 

k 

ki^hk-Q 

1+7.961 

1+7.918 

1+7 .  801 

5 

h-k^^os 

1+5,181 

1+5.11+1 

i+l+.958 

6 

i|l.I^89 

1+2.365 

1+2. 327 

1+2.096 

0 

1U;.629 

i+l+.712 

1|1|.682 

l+l+.668'^ 

1 

1+2.665 

1+2.790 

i|2.76l 

1+2.71+2 

2 

1+0.731+ 

1+0.905 

1+0.877 

1+0.838 

HD        3 

38.833 

39.01+7 

39,020 

38.998 

1; 

36,923 

37.251 

37.226 

37.11+0 

5 

35.026 

35.1+21 

35.397 

6 

33.087 

33.620 

33.597 

0 

29.890 

29.91+3 

29.930 

29.910® 

1 

28.816 

28.891+ 

28.881 

28.81+8 

2 

27.759 

27.857 

27.81+5 

°2        3 

26,716 

26.838 

26.826 

i; 

25.681 

25.81+7 

25.836 

5 

2i+,627 

21^.81+0 

2i+.  829 

6 

23.896 

23.885 

Calculated  with  McLean,    Weiss,    and  Yoshlmlne's   inter- 
nuclear  potential    [1+7]. 

Calculated  with  Kolos    and  Roothaan's   internuclear 
potential    [[|.8J. 

Herzberg   and   Howe    [51]. 

Durie    and  Herzberg    [50]. 

®Stoicheff    [1+9], 
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TABLE  10.   CALCULATED  AND  OBSERVED  ROTATIONAL  CONSTANTS 

D^   AND   H    FOR   H    HD,   AND  D^       (IN  CM""'-). 


MWY" 


KR^ 


Exp. 


V 

D 

V 

H^xlO^ 

D 

V 

H^  X  10^ 

D 

7 

H^x  10^ 

0 

.0i|56 

2.i| 

,014.56 

3.6 

.01+59° 

.  01+59 '^ 

S,2^ 

1 

♦  Oi^i^ 

5.2 

.0104-3 

5.1+ 

.01+32 

.01+37 

2 

,01+30 

3.0 

,014.27 

1+.8 

.01+27 

.01+27 

H2     3 

.014.16 

l+,0 

,01+11+ 

5.7 

.01+12 

.01+09 

3.75^ 

k 

.01+06 

.01+01+ 

.0397 

3.50 

5 

.0399 

.0389 

.0385 

3.1; 

6 

.0385 

.0377 

3.3 

0 

,0259 

3.5 

.0257 

1.9 

.0259° 

.0263"^ 

2.2° 

1 

.0250 

lol 

.0251 

2.6 

.0255 

.02514- 

2 

.02l|l; 

l.ll 

,021+3 

2,7 

.02)11 

HD      3 

.0237 

2,0 

.0237 

2.3 

.021+1+ 

k 

00232 

,021+2 

.0231 

5 

.0227 

,0221+ 

6 

.021+2 

0 

,0111; 

.0111+ 

-.1 

.0113° 

.36° 

1 

.0112 

.0113 

,9 

,0107 

2 

.0110 

.0110 

1.6 

^2     3 

.0107 

,0108 

1.2 

i; 

.0105 

.0106 

5 

.0102 

.0110 

6 

.0109 

Calculated  with  McLean,  Weiss,  and  Yoshimine ' s  inter- 
nuclear  potential  [I+7J. 

■J- 

Calculated  with  Kolos  and  Roothaan's  internuclear 
potential  [1+8J. 

°Stoicheff  [1+9]. 
^Herzberg  and  Howe  [5l]» 
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The  discrepancy  of  -.012%  In  the  equilibrium  values, 

R  =  1.1^0100  and  R  =  I.I4.OO83,   of  the  MWY  and  KR 
e  6 

potentials,  respectively,  would  correspond  to  a  discrepancy 
of  .006%  in  B   for  the  rigid  rotator.   Consequently, 
the  difference  of  approximately  <.2/o  in  the  calculated  B^ 
values  must  be  due  to  the  poorer  anharmonicity  of  the  MWY 
potential. 

The  experimental  and  calculated  data  for  the   v  =  0 
vibrational  level  of  the  heaviest  isotope,   Dp,   should 
provide  the  best  comparison  of  the  theoretical  and  experi- 
mental  R    since  this  is  the  case  where   R   is  nearest  to 
e 

R   durine;  the  vibrational  motion  and  where  Van  Vleck's 
e       ^ 

corrections  are  smallest  and  most  accurate.   In  Tables  8 
and  9  it  is  seen  that  the  calculated  S  (j)   and   B 
values  are  approximately   .06%   higher  than  the  experi- 
mental ones.   This  implies  that,  if  one  accepts  the 

Van  Vleck  correction,  the  calculated   R    of  the  KR 

'  e 

potential  is   »03%>   lower  than  the  experimental  evidence 
indicates  it  should  be.   However,  the  figure  of   .O6/0   is 
so  close  to  the  Van  Vleck  correction  of   .OI4J4./0   to   B 
and  to  the  numerical  and  experimental  errors  that  this  can 
only  be  accepted  as  an  order  of  magnitude  estimate. 

The  calculation  of  the  anharmonicity  corrections 
defined  in  (V  I.2I4.)  is  described  in  Table  11.   The  value 
of   Yj   which  is  used  to  calculate  <^",   was  obtained  by 
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TABLE  11.   CALCULATION  OF  ANHARMONICITY  CORRECTIOr 


^2 

HD 

^2 

D^^ 

38  285.59 

38  285.59 

38  285.59 

^->-> 

38  279.06 

38  281,56 

38  28[|.23 

A. 

6,53 

k'03 

1,36 

e 

60,913 

14-5.699 

30.[;80 

Be 

60.905 

[|5. 690 

30,1|70 

Ya 

-„00013 

-,00020 

-,00033 

-: 

iUl.10.99 

3820.51 

3120.21; 

u.- 

)|)|12.67 

3823^51^ 

3125  ..oil 

\ 

-,00038 

-,00079 

-,ooi5i| 

S-r 


The    "''«■"    denotes    quantities   derived   from   the 

equilibrium    constants      D   =    .17i|-  144-2   a,u., 

Rg   =  l,lj.0083   acU»,    and     y  =    O7083   a,u.    obtained 

from  the   Kolos    and   Roothaan  potential    [1|8J. 

The    "^H^"    denotes   quantities   obtained  by   an  exact 

fit   of   the   energy  formula  to   the    calculated 

V  =  0,1,2,3,      J  -   0,1, ..,,14.     term  values. 


numerical  differentiation  of  the  slope  of  the  internuclear 
potential  as  given  by  the  Virial  Theorem, 

The  calculated  molecular  constants  obtained  with  the 
KR  potential  are  siimmarized  in  Table  12,   The  calculated 
quantities  in  parentheses  were  obtained  by  applying  the 
Van  Vleck  corrections  (2,1).   The  experimental  quantities 
in  parentheses  are  of  doubtful  accuracy  due  to  a  lack  of 
experimental  evidence.   It  must  be  remembered  that  all 
calculated  quantities  with  the  subscript  "e"  are,  for  the 
reasons  given  above,  subject  to  some  uncertainty  due  to 
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TABLE  12.   VIBRATIONAL-ROTATIONAL  CONSTANTS  OP 
H„,   HD,   AND   D„   (IN  CM"-"- )  . 


HE 

°2 

Gale.             Exp.^'^ 

Gale. 

Exp.^'^ 

Gale, 

Exp,^'^ 

B 
o 

59.398 
(59.3i;5) 

59.339 

,14^.712 
([|.i+.682) 

)|)|,669 

29.91+3 
(29,930) 

29,911 

B 
e 

60.905 
(60,851.) 

60.861^ 

1+5.690 
(1+5.659) 

1+5.663 

30,1+70 
(30.1+57) 

(30,1+1+2) 

a 
e 

3.037 

3.076  1+ 

1.970 

2,003 

1,051+ 

Ye 

.01^.988 

.060  1 

.026   8 

.039   7 

,000  6 

5 
e 

-.003  5 

-.001  8 

-.003  h 

+  .001 

R 
e 

.7i|l  66A 

,71+1  58A 

.71+1  58A 

.71^1  51A 

.7l;l  57A 

.71+1  65A 

r: 

o7!|l    28A 

.714-1   16A 

,714-1   28A 

.71+1  19A 

,7'n    28A 

.7l;l  1+3A 

D 
e 

.Oi;5  92 

.01+6  57 

.025  9 

.026   7 

,011  1+ 

(.011  61+) 

Pe 

-.000  [j.6 

-.001   62 

-.000   1 

-.000  7 

.000  05 

H 
e 

l|5xl0"^ 

1.5x10"^ 

2.2x10"^ 

-,5x10"^ 

-T 
oo 

36  103.5 
(36  109.0) 

36  113.0+.3 

36  393.2 
(36  397.8) 

36  399,9+1, 

36  737.7 
(36   71+0,5) 

36  7l;3«6+o5 

D® 

38  279.1 
(38  281^0 6) 

38  292.31-5 

38  281.6 
(38  286,2) 

38  290,3+1. 

38  281+.2 
(38  287.0) 

38  290,8+,7 

VG(1) 

i;l69.29 
(1+166.33) 

1+161.11+ 

3638,98 
(3637.09) 

3632,15 

2999,02 
(2997.97) 

^e 

i+l+12o67 
( 14409. 5i|) 

1+1^00.39 

3823.51+ 
(3821.55) 

3812.29 

3125,01+ 
(3123,95) 

'^e^e 

123»71 

120,82 

93.935 

90,908 

61+,  1+68 

^^^e 

lc2[+3 

.721^ 

1.017 

.501+ 

,898 

Herzberg  and  Howe  [5lJ,  with  the  exception  of  -T    and  D. 

Durie  and  Herzberg  [50j,  with  the  exception  of  -T   ,D,andH 
The  latter  is  given  by  Stoicheff  [l+9j.  °°        ® 

°Stoicheff  [1+9  J,  with  the  exception  of  -T    and  D. 

d  °° 

The    experimental    -T        and   D  values   are    those   reported  by 

Herzberg    and  Monfils    [52 J, 

The    calculated  D-values    given  were    obtained  by   fitting   the 
energy  formula   to    the    calculated   term  values.      The   value    of 
D  from  Kolos    and   Roothaan's   potential    ciirve    is 
D  =  38  285,6   cm"^. 
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the  ambiguity  in  their  def mititi ons .   In  the  "Calc." 

columns,   R"  was  obtained  directly  I'rom  the  KR  potential 
'    e 

and   R   was  obtained  from  the  calculated  B    after  the 

Van  Vleck  correction.   The  values  of  R^   and  R"   in  the 

"Exp."  coliOTins  are  those  given  by  the  authors  of  the 

experimental  data  and  are  to  be  compared  with  the  calcu- 

lated   R   and   R  „ 
e        e 

It  may  be  observed  in  Table  12  that  the  corrected 
theoretical  energies   -  T^^   of  E^,    HD,  and  D^   are 
I;,  2,  and  3  cm.~   lower  than  the  corresponding  experi- 
mental values.   However,  the  potential  used  here  was  for 
the  [;0-term  electronic  wave  function.   A  50-term  wave 
function  was  also  calculated  by  KR  for   R  near   R   and 
gave  a  value  of   D  which  was  1,5  cm.    higher.   Adding 

this  figure  to   -T    reduces  the  discrepancy  with 
^  oo 

experiment  almost  to  numerical  and  experimental  error. 
As  Herzberg  and  Monfils  [S2]    point  out,  FrBman's  [531 
estimates  of  relativistic  effects  indicate  that  one  may 
have  to  add  between  3  and  23  cm,    to  these  theoretical 
values  while  ele ctrodynamic  (Lamb  shift)  corrections 
would  require  the  subtraction  of  some  small  amount. 

These  results  indicate  that,  in  addition  to  giving 
an  accurate  binding  energy,  the  KR  potential  gives  a 
correct  prediction  of  the  vibrational-rotati onal  levels 
and  consequently  has  the  correct  shape  for  a  wide  range 


-  145 


of  R-values.  This  gives  one  confidence  in  the  accuracy  of 
the  nuclear  wave  functions  obtained  and  used  for  calculat- 
ing the  expectation  values  reported  in  the  next  section. 

3.   Expectation  Values  ^R  )  ^  . 

The  expectation  values  of  powers  of  internuclear 
distance  were  obtained  from  the  solutions  of  the  nuclear 
wave  equation  with  the  KR  internuclear  potential.   The 
results  are  given  in  Tables  13-15  in  a.u.   It  is  expected 
that  the  use  of  these  results  will  yield  some  improvement 
in  the  derivation  of  important  physical  constants  from 
radiofrequency  spectra. 

One  simple  and  direct  application  can  be  made  by 
referring  to  Ramsey's  [Sh-]   measvirement  of  the  nuclear 
spin-spin  interaction  energy 

(3.1)  d  =  ^,2(r-3>3^, 

of   Hp.   Taking,  for  the  proton  magnetic  moment, 

[i-p  =  1.52102x10  -^     Bohr  magnetons,  Ramsey's  experimental 

findings  and  the  present  theoretical  results  give 

<R''^^q  1  ^  Q'^SS   873  and  0.356  099  a.u. 
respectively. 
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TABLE  13.   EXPECTATION  VALUES  <(R^/^^j   for   H^   (IN  A.U.  ) 


V 

/ 

J 

0 

1 

2 

0 

-3 

.357  756 

,356  099 

,352   821; 

-2 

.14.96  928 

.1+95  1+014. 

,1;92   386 

-1 

.700  061; 

o698  999 

.696   882 

1 

l.)|)|7  989 

1,1+50  165 

1,1+51+  505 

2 

2. 1214.  832 

2,131  178 

2.11+3   861+ 

1 

-3 

.3kS  068 

,31+3  1+19 

,31+0  160 

-2 

.i|72  2i;7 

.1+70  772 

,1+67  851 

-1 

.673  263 

,672  231+ 

,670  189 

1 

1.5i|4  370 

1.51+6  639 

1,551  166 

2 

2.i;70   814.1; 

2, 1^77   980 

2,14.92  21|8 

2 

-3 

.331   786 

.330  151+ 

.326   932 

-2 

.)lll8   11;0 

.)ll|6  716 

.)|)|3  896 

-1 

.61^7  157 

.6I4.6  163 

,6l4.)|  189 

1 

1.6144  722 

1.61+7  092 

1.651  822 

2 

2.852   156 

2.860  175 

2.876   211 

3 

-3 

.317   837 

.316  221+ 

.313   OI4.O 

-2 

.I1.2I;  1^22 

.1+23  01+3 

.I4.20  313 

-1 

.621  558 

.620  593 

,618  679 

1 

l,7l|9   81^8 

1.752   31+3 

1.757  321 

2 

3. 271;  379 

3.283  1+36 

3,301  51+6 
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TABLE   11;.      EXPECTATION    VALUES 


<R^>, 


vJ 


FOR      HD      (IN  A.Uo ) 


V 

/ 

J 

0 

1 

2 

0 

-3 

.358  558 

.357  313 

.351;  8i;7 

-2 

,1^98  620 

.1|97  1+71+ 

.1+95  197 

-1 

,701   908 

.701  106 

,699  511 

1 

1,^1+1  591+ 

1.1443   222 

1.10+6  1+72 

2 

2,102  14.514- 

2.107   177 

2,116   620 

1 

-3 

.31+7   653 

.31^.6  14.13 

o31+3  955 

-2 

.1+77   186 

.I4.76  071 

.1+73  857 

-1 

0678   618 

.677  839 

,676  291 

1 

1,521;  611| 

1.526  302 

1,529  672 

2 

2o398   180 

2.I1O3    1;09 

2  c  1+13  866 

2 

-3 

.336  267 

.335    039 

,332  605 

-2 

.I4.56   168 

.1;55  088 

^h^2  91+1+ 

-1 

,655  850 

.655  097 

.653  597 

1 

1,610  51;9 

1,612   298 

I.615  790 

2 

2,719  752 

2.725  527 

2.737  076 

3 

-3 

.321;  350 

,323  129 

.320   713 

-2 

.1;35  1^0 

.1;31+  387 

.1+32  296 

-1 

.633  1;72 

.632  736 

.631  273 

1 

1,699  952 

1.701  782 

1.705  1+31+ 

2 

3,070  828 

3,077  253 

3,090  099 
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TABLE   15.      EXPECTATION    VALUES     (R'^^^t      FOR     ^p      (IN  A.U.  ) 


V 

X 

J 

0 

1 

2 

0 

-3 

.359  506 

.358  675 

.357   021+ 

-2 

.500  625 

.i|99  858 

.1+98   331 

-1 

,70^  091^ 

.703  557 

.702  1+87 

1 

l.i|3l+  01^8 

i.i;35  131 

1.1+37   291+ 

2 

2.076  185 

2.079  306 

2,085  51+6 

1 

-3 

.350  705 

.3i;9  875 

.31+8   227 

-2 

.I483   082 

.14-82  332 

.1+80   838 

-1 

.685  010 

.681;  1+86 

.683   l|'|2 

1 

1.501  385 

1.502  501 

1,501+  729 

2 

2.313  853 

2.317   21+7 

2.321+  033 

2 

-3 

.31^-1  i|75 

.31+0  653 

.339   021 

-2 

.1;65  737 

.J4.65  007 

.1+63  555 

-1 

.666  228 

.665  719 

.661+  705 

1 

1.570  716 

1.571  860 

1.571+  11+6 

2 

2,568  558 

2.572   228 

2.579  561+ 

3 

-3 

,33>2  006 

.331   183 

.329  51+7 

-2 

.)|)|8  711 

.)|)|7  991 

.)|)|6  560 

-1 

.6I4.7  805 

.61+7  302 

.61+6  302 

1 

i.6;i2  079 

1.61+3  273 

1.61+5  657 

2 

2,8i;l   1;80 

2.81+5  507 

2.853  560 
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Ll.   Expectation  Values  \   >   r  ^~  nl  P  (cos  9  .)/  ^. 
^       — ^ \  '  —->      ai       n      ai  /vJ 

The  expectation  values  considered  here  give  the 

average  values,  over  electronic  and  nuclear  coordinates, 

of  the  operators 

(1,.2)  5:  r;^  cos  9^^ 

(i+.3)  YZ  ^ll   (3  cos^  9^^-  1) 

where  the  subscript   ai   refers  to  the  spherical  polar 
coordinates  of  electron  i  referred  to  nucleus  A.   The 
calculated  values  listed  in  Table  16  were  obtained  with 
the  MWY  electronic  wave  function  and  internuclear  potential 
which,  it  should  be  remembered,  is  less  accurate  than  the 
KR  potential  used  in  the  previous  sections.   Expectation 
values  of  (l+.l)  and  (I|..3)  give  the  nuclear  screening 
parameter  and  electric  field  gradient  required  in  the 
analysis  of  radiofrequency  spectra.   The  average,  over 
electronic  coordinates,  of  {l\..2)    is  the  electric  field  at 
nucleus  A  which,  by  the  Hellman-Feyninan  Theorem,  is 
approximately  equal  to  the  force  on  A  due  to  the  electrons 
at  the  distance   R   at  which  (L|..2)  is  calculated.   The 
average  of  this  quantity,  over  R      is  then  simply  an  estimate 
of  the  expectation  value  of  the  electric  field  at  A, 
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2  ,  . 

TABLE  16.      EXPECTATION    VALUES   <^  ^^  r^J        ni    P^(cos    ©^1  VvJ 


2_ 
i~l 


FOR     H^,      HD,      AND     D^      (INA.U.). 


V 

n 

0 

J 

1 

2 

0 

0 
1 
2 

1,799   119 
,[|9i|  280 
.359  572 

1,797   829 
.1^93  515 
.358  71fO 

1.795  265 

.  .1;91  993 

.357  089 

1 

TT 

0 
1 
2 

lo760   313 
ol+69   717 
.338  275 

1.759   086 
.1;68  983 
.337  1;89 

1,756   61;9 
.1;67  523 
.335  928 

II2        2 

0 
1 
2 

1,722  571 
^hk^  560 
o317  37IJ. 

1,721  387 
.316   626 

1,719   032 
.1+1;3  1;19 
.315  lUO 

3~^ 

0 
1 
2 

1.681;  136 
,1+20  307 
.296  563 

1.682  957 
.1;19  576 
.295  831 

1,680  611 
,I;18   119 
.291;  377 

0 

0 

1 
2 

1.801  817 
.[;95  978 
.361   061 

1.800  8);6 
.i;95  1^02 
.360  1+31; 

1,798   912 
.1;9);  255 
,359   187 

1 

0 
1 
2 

1.768  01^1 

^klk  633 
.31^2  563 

1.767   112 
.1^714-  078 
.3lil  967 

1.765   263 
.1;72   972 
.31^0   781 

XliJ           p      - 

0 
1 
2 

1.735  229 
.1^53  708 
>32k  321 

I.73I;  337 
.lj-53  172 
.323  751; 

1,732  560 
.1;52   101; 
,322  625 

3 

0 
1 
2 

1,702  I9i| 
.[4.32  296 
.306   291 

1,701  303 
.1+31  71;8 
.305  731f 

1,699  52B 
.1;30  655 
.301;  626 

0 

0 
1 
2 

1,805  019 
,1;97  989 
.362   820 

1,801;  368 
.1;97    601; 

.362  a.00 

I.803   069 
,1;V6  831; 
.361   562 

1 

0 
1 
2 

1.777   287 
,[|80  505 
.3^7   6914. 

1.776  659 
.1;80   131 
.31^7    290 

1.775  1;07 
.1;79  381; 
.31^6  1;86 

^2        2 

0 
1 
2 

1.750  205 
,1|63  257 
.332  601 

1.7l;9  603 
.1;62   898 
.332  217 

1.71;8  1;01; 
.1;62   180 
.331  l;5o 

3 

0 

1 
2 

1.723  1;73 
,)|!|6   210 
.317   821 

1.722   880 
.1015   81|9 
.317  1U4-5 

1.721   696 

.hhS  131 
.316   696 
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The  measured  value  ISh-^  of  the  quadrupole  coupling 
constant  in  the  v  =  0,  J  =  1  level  of  the  electronic 
ground  state  of   Dp   is 

ik'h-)  eQq   =   (1. 1^907  +  .0015)  x  lO"^^  ergs 

where   Q   is  the  deuteron  quadrupole  moment  and 

q   =   <2R-^>o  1  -  Q,, 

(1^.5)  2 

Q,,  =  <^''li   (3  cos2  e^,-i)>o^,  . 

Newell  [55J  calculated   q  =  ,3S2     which  is  close  to  the 
value   q  =  '35^      obtained  from  the  data  in  Tables  15  and 
16.   Newell 's  value  of   q   has  been  used  with  (14-. I4.)  to 
obtain  the  accepted  value  of  Q.   Recently,  Auffray  [56 j 
obtained   q  =  .339  with  a  wave  function  consisting  of 
10  of  the  L|.0  terms  used  by  KR.   He  performed  a  more 
thorough  optimization  of  the  scale  factor  than  KR  did  and 
obtained  approximately  the  same  electronic  energies.   It 

was  shown  in  Section  VI  3  that  the  LCAO-type  wave  function 

+ 
for   H„   gave  values  of   Q,    which  were  too  low.   One  can 
d  zz 

infer  that  a  similar  deficiency  will  occur  in  the  MV\iT  wave 
function  and  that  Auffray 's  result  at  least  gives  a  change 
in  the  right  direction. 
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